R Brenner

Numerical computation of the response of piezoelectric composites using Fourier transform

Keywords: numbers: 46.25.Hf, 47.11.Kb composite, piezoelectric, Fourier transform, periodic homogenization

An efficient numerical scheme is presented to compute the response of piezoelectric composite materials with arbitrary complex microstructures. It makes use of Fast Fourier Transform (FFT) to solve iteratively the coupled periodic Lippmann-Schwinger equations for an heterogeneous electroelastic medium. The method is assessed in the case of a two-phase composite by comparison with analytic solutions and finite-element results taken from the literature.

I. INTRODUCTION

This work aims at describing an efficient numerical scheme to determine the overall behaviour and the local fields within an heterogeneous medium presenting a coupled constitutive behaviour. The method is an extension of a computational approach initially proposed for elastic-type composites 1 . Up to now, numerical homogenization in the framework of coupled behaviours (e.g piezoelectric, piezomagnetic, thermoelastic etc.) has been essentially performed with the Finite-Element method (FEM) 2-5 which requires tedious meshing as soon as the microstructure becomes complex. The alternative method presented in this article has been previously used for uncoupled (linear and nonlinear) mechanical behaviours with various microstructures 1, [6][7][8] as well as other uncoupled physical properties 9,10 . It makes use of Fast Fourier Transform (FFT) and permits to use images of the microstructure for the computation.

Section II is devoted to the mathematical formulation of the problem to be solved. It is shown that the local fields are solutions of a set of coupled Lippmann-Schwinger (LS) equations. The iterative scheme which allows to solve the LS equations is detailed in section III. Results are presented for a two-phase composite material with laminate or fibrous microstructure in section IV.

II. INTEGRAL EQUATIONS FOR HETEROGENEOUS ELECTROELASTICITY

We consider a Representative Volume Element (RVE) Ω of a heterogeneous electroelastic material under isothermal conditions. The local properties of the material are thus defined by the elastic C, the dielectric (permittivity) γ and the piezoelectric e moduli tensors.

Relevant field variables to describe the behaviour of the medium are thus the elastic strain ε, the stress σ, the electric induction D and the electric field E. The elastic strain and the electric field are respectively derivable from the displacement field u and the electric potential φ. Tensor rank and units of the different variables of the problem are summarized in Table I.

To solve the problem, periodic boundary conditions are chosen and the RVE is considered as a unit cell. Using periodicity conditions stated by Suquet 11 , the fields u and φ on the

σ 2 N m -2 Elastic strain ε 2 - Elastic displacement u 1 m Electric induction D 1 C m -2 Electric field E 1 V m -1 Electric potential φ 0 V Piezoelectric moduli e 3 C m -2 Permittivity moduli γ 2 F m -1 Elastic moduli C 4 N m -2
boundary ∂Ω of the unit cell are decomposed into

u(x) = ε.x + u ⋆ (x) and φ(x) = -E.x + φ ⋆ (x) (1)
with ε and E prescribed average values of the strain and electric fields over Ω. u ⋆ and φ ⋆ are periodic (#) fluctuating terms. Besides, the mechanical and electric static equilibriums between neighboring cells imply anti-periodic (-#) conditions for the traction vector σ.n and the electric charge D.n on ∂Ω (n is the outer unit normal on the boundary).

The problem to be solved can thus be expressed as: find u(x) and φ(x), ∀x ∈ Ω, such

that                    Potential fields: ε = 1 2 (grad u + T grad u), E = -grad φ Local laws: σ(x) = C(x) : ε(x) -T e(x).E(x) D(x) = e(x) : ε(x) + γ(x).E(x)
Static equilibrium:

div σ = 0, div D = 0
Boundary conditions: u ⋆ , φ ⋆ # and σ.n, D.n -#, ∀x ∈ ∂Ω.

(2)

Note that the transpose of e is defined by T e ijn = e nij . The notation introduced by Barnett and Lothe 12 to treat electroelastic problems is unnecessary for our study. Simple tensorial relations are kept to describe the problem.

By introducing a homogeneous reference medium with elastic C 0 , permittivity γ 0 and piezoelectric e 0 moduli tensors, the heterogeneous electroelastic problem is transformed into a homogeneous electroelastic problem with heterogeneous fields of fictitious body forces (f ′ ) and free electric charges (q ′ ) 13,14 . Static equilibrium relations thus read div(C 0 : grad u + T e 0 .gradφ) + f ′ = 0 and div(e 0 : grad uγ 0 .gradφ) + q ′ = 0 (3) with f ′ = div(δC : grad u + T δe.gradφ), q ′ = div(δe : grad u -δγ.gradφ).

The notation δa means a(x)a 0 . The system of differential equations (3) can be solved using standard Green's functions method. In linear piezoelectric media, the electroelastic Green's functions describe the elastic displacement u and electric potential φ at point x due to a unit point force or electric charge at point x ′ . The solution to this coupled-field problem was derived by Deeg 15 . It has been used to obtain a generalization of Eshelby's inclusion solution in elasticity 16,17 and to extend different micromechanical models to coupled electroelasticity [17][18][19] . By making use of the superposition principle, Deeg's result allows to get the solution of the periodic homogeneous electroelastic problem (3) with fluctuating fields f ′ (x) and q ′ (x). It is expressible as coupled periodic Lippmann-Schwinger equations

ε(x) = ε -Γ 0 * τ (x) -T Υ 0 * P (x), ∀x ∈ Ω, E(x) = E + Υ 0 * τ (x) + ∆ 0 * P (x), ∀x ∈ Ω, (5)
which are linear integral equations for the fields ε and E (the notation * indicates a convolution product). τ and P are polarization fields from which f ′ and q ′ are respectively derivable. They read τ = δC : ε -T δe.E and P = δe : ε + δγ.E.

It can be remarked that only three Green's operators Γ 0 , ∆ 0 and Υ 0 have to be introduced because of the symmetry of the piezoelectric effect (i.e the value of the strain ε at point x due to a unit negative electric induction at point x ′ is equal to the value of the electric field E at point x due to a unit stress at point x ′ ). The coupled periodic Lippmann-Schwinger equations ( 5) read in Fourier space

ε(ξ) = -Γ0 (ξ) : τ (ξ) -T Υ0 (ξ). P (ξ), ∀ξ = 0, ε(0) = ε, Ê(ξ) = Υ0 (ξ) : τ (ξ) + ∆0 (ξ). P (ξ), ∀ξ = 0, Ê(0) = E. (7)
The operators Γ 0 , ∆ 0 and Υ 0 are explicitly known in Fourier space for whatever anisotropy of the reference piezoelectric medium 16,19 . They are given by

Γ0 ijkl (ξ) = 1 4 K -1 il ξ k ξ j + K -1 jl ξ k ξ i + K -1 ik ξ l ξ j + K -1 jk ξ i ξ l , Υ0 nij (ξ) = 1 2 K -1 4j ξ i ξ n + K -1 4i ξ j ξ n , ∆0 ij (ξ) =K -1 44 ξ i ξ j , (8)
with K a 4 × 4 matrix whose components read

         K ij =C 0 pijq ξ p ξ q , ∀(i, j) ∈ [1, 3] 2 , K 4i =K i4 = e 0 piq ξ p ξ q , ∀i ∈ [1, 3], K 44 = -γ 0 pq ξ p ξ q . (9)
When the reference medium does not present piezoelectric coupling (i.e e 0 = 0), Υ0 = 0 while the operators Γ0 and ∆0 reduce, respectively, to the Green operators for elasticity and electrostatics. It is stressed out that the presented formalism can be applied to other coupled-field problems which have the same structure (i.e constitutive laws, divergence and gradient equations). In particular, it can be used to tackle magnetoelectroelastic problems.

III. FOURIER TRANSFORM BASED METHOD

A. Iterative scheme

An efficient numerical scheme based on Fourier transform has been previously proposed to solve Lippmann-Schwinger type equations for elastic 1 , conduction 9 or dielectric 10 properties.

This method uses the fact that the polarization fields are defined locally in real space whereas the Green's operators act locally in Fourier space where they are explicitly known.

The iterative scheme proposed by Moulinec and Suquet for elastic problems 1,20 can be extended to the context of composite materials with coupled physical properties. In the case of piezoelectric materials, the direct Fourier transform of an initial guess of the polarization fields τ and P (6), defined on a regular grid, allows to obtain the strain and electric fields in Fourier space. By applying inverse Fourier transform to these latter, the stress and electric induction fields can be computed in real space using the constitutive equations and lead to a new estimate of the polarization fields. The solution fields ε and E are then obtained by successive iterations until convergence criterion is reached. Details on the rate of convergence of the method can be found in Michel et al. 21 .

An attractive feature of this approach is that the pixel discretization of the image of the microstructure can be used to assign the local properties of the material at the nodes of a regular grid. Compared to the FEM, there is thus no need to mesh the microstructure. This allows to consider easily materials with complex microstructures 1,6 . Besides, by using Fast Fourier Transform (FFT) algorithms, the computational cost scales almost linearly with the size of the problem (i.e number of pixels in 2 dimensions and voxels in 3 dimensions).

B. Algorithm

Periodic microstructures are considered. Consequently, the microstructure of the material is obtained by repeating periodically an elementary pattern (i.e unit-cell). First, this unitcell is discretized in a regular grid and the local properties C, e and γ are assigned at each point x g of the grid. Then, the reference medium has to be defined to compute the different Green operators (Γ 0 , Υ 0 , ∆ 0 ) in Fourier space at each frequency ξ g . The convergence of the method is strongly sensitive to the choice of the reference medium. Based on previous results in the mechanical context 1,21 , we chose C 0 = C , e 0 = e and γ 0 = γ .

Following Moulinec and Suquet 1 , the iterative algorithm then reads as follows:

Initialization: ε(x g ) = ε, E(x g ) = E, ∀x g Iterate: 1. σ(x g ) = C(x g ) : ε(x g ) -T e(x g ).E(x g ) D(x g ) = e(x g ) : ε(x g ) + γ(x g ).E(x g ) if convergence reached -→ stop 2. τ (x g ) = σ(x g ) -C 0 : ε(x g ) + T e 0 .E(x g ) P (x g ) = D(x g ) -e 0 : ε(x g ) -γ 0 .E(x g ) 3. τ (x g ), P (x g ) F -→ τ (ξ g ), P (ξ g ) 4. ε(ξ g ) = -Γ0 (ξ g ) : τ (ξ g ) -T Υ0 (ξ g ). P (ξ g ), ∀ξ g = 0, ε(0) = ε Ê(ξ g ) = Υ0 (ξ g ) : τ (ξ g ) + ∆0 (ξ g ). P (ξ g ), ∀ξ g = 0, Ê(0) = E 5. ε(x g ), E(x g ) F -1 ←-ε(ξ g ), Ê(ξ g )
Convergence is reached when the stress σ and electric induction D fields are in equilibrium within a certain tolerance error. The convergence criterion is defined as

max ξ g . σ(ξ g ) 2 1 2 σ(0) , |ξ g . D(ξ g )| 1 2 D(0) ≤ η. (10)

IV. ILLUSTRATIVE APPLICATIONS

To assess the proposed numerical scheme, we consider a binary composite made of an isotropic polymer phase (epoxy) and a transversely isotropic piezoelectric ceramics (lead zirconium titanate -PZT) since numerous results exist in the literature. The constitutive law of materials which are transversely isotropic about the e 3 direction can be written, using Kelvin's notation (see Appendix A), as

                      σ1 σ2 σ3 σ4 σ5 σ6 D 1 D 2 D 3                       =                       C11 C12 C13 0 0 0 0 0 -ȇ 31 C12 C11 C13 0 0 0 0 0 -ȇ 31 C13 C13 C33 0 0 0 0 0 -ȇ 33
0 0 0 C44 0 0 0 -ȇ 15 0 0 0 0 0 C44 0 -ȇ 15 0 0 0 0 0 0 0 C66 0 0 0 0 0 0 0 ȇ15 0 γ 11 0 0

0 0 0 ȇ15 0 0 0 γ 11 0 ȇ31 ȇ31 ȇ33 0 0 0 0 0 γ 33                       .                       ε1 ε2 ε3 ε4 ε5 ε6 E 1 E 2 E 3                      
The isotropic constitutive law is obtained with the following additional relations

C33 = C11 , C12 = C13 , C66 = C44 , ȇ = 0, γ 33 = γ 11 . (11)
Material data for PZT and epoxy are given in Table II. The FFT-based numerical scheme is now used to compute the response of the binary composite with two different arrangements of the constituents. The simplest microstructural case is the laminate composite which is made of succcessive layers of the constituents with arbitrary thickness. Indeed, it is well known that there is an exact solution to the problem and that the fields are constant within each layer. The solution has been given by Benveniste and Dvorak 22 in the specific case where the axes of transverse isotropy of each layer are parallel to each other and perpendicular to the lamination direction. In the sequel, e 3 is the axis of transverse isotropy and e 2 is normal to the layers.

The applied mechanical loading is shear parallel to the layers (ε 23 imposed) and the electric loading is normal to the layers (E 2 imposed). The nonvanishing components of the strain ε and electric field E within layer (r) read 22 ε (r)

23 = A (r) ε ε 23 + B (r) ε E 2 , E (r) 2 = A (r) E ε 23 + B (r) E E 2 . ( 12
)
For completeness, the expressions of the influence factors

A (r) ε , B (r) ε , A (r) E and B (r) E are given in Appendix B.
The results of the numerical scheme are presented in Figure 1 for the local strain and electric fields within a binary laminate composite with c 1 = c 2 = 0.5. The unit-cell has been discretized with a regular grid of 32 × 32 pixels and the convergence error η is equal to 10 -6 .

Only two iterations are necessary to reach convergence in this case. The FFT computational scheme leads to uniform fields within each layer and the analytical solution ( 12) is obtained numerically. It is worth noting that convergence is obtained even if the composite present an infinite piezoelectric contrast because of the finite contrast on the elastic and dielectric properties.

B. Fiber composite

We further consider a composite material made of PZT parallel circular fibers embedded in an isotropic epoxy matrix. In the cross section of the material, the fibers are arranged on a periodic square array. A representative unit-cell is thus a square with a fiber at the center and the overall behaviour exhibits quadratic symmetry.

To compute the effective moduli of the composite, the unit-cell has been discretized with a regular grid of 512 × 512 pixels and we applied successively six different loadings (ε 11 , ε 33 , ε 23 , ε 12 , E 1 and E 3 ). The results obtained with the FFT-based numerical scheme are compared with previously published FEM studies 2,3 in Table III. For this comparison, note that we converted the overall moduli obtained by FEM which were given at constant electric induction D whereas we used, in the present work, constitutive equations at constant electric field E. An excellent agreement is obtained with the FEM results for all effective moduli and the computational method requires, in average, only a few iterations to converge with an error of 10 -6 on the fields equilibrium. Besides, it has been checked that our results fulfill, with a relative error of a few percents, the universal relations for the overall moduli of fiber composite materials with quadratic (tetragonal) symmetry 22 . For illustrative purpose, the distribution of the local fields E 1 and ε 13 for a macroscopic loading E 1 are plotted in Figure 2. Almost uniform fields are predicted within the fibers where, as expected, a low electrical field takes place. Significant longitudinal shear strain is predicted in the matrix at the fiber-matrix interface normal to the e 1 axis. By following the approach initially developed for uncoupled mechanical constitutive behaviour by Moulinec and Suquet 1 , we have presented an original computational scheme to predict the local and overall responses of heterogeneous electroelastic materials. This scheme, which makes use of fast Fourier transform algorithms and allows to use the image of the material as direct input, is a viable alternative to finite-element analysis to tackle the response of piezoelectric composites. In this article, the method has been successfully compared to analytical solutions and finite-element results for simple binary composites.

However, the method is particularly well suited to consider heterogeneous materials with complex microstructures like polycrystalline ceramics. Such microstructural computations have been previously performed in the uncoupled mechanical context 6,8 . Eventually, it is noted that this computational procedure can be used for nonlinear heterogeneous materials with coupled constitutive equations and is not restricted to the coupled piezoelectric behaviour. It is relevant for other "multiphysics" problems which adopt a similar mathematical description (e.g magnetoelectroelasticity).

Appendix A: Kelvin notation for symmetric tensors

The Kelvin notation 23 , initially introduced in elasticity, can be advantageously used to write the electroelastic constitutive equations as

   σα = Cαβ εβ -T ȇαj E j D i = ȇiβ εβ + γ ij E j (

FIG. 1 :

 1 FIG. 1: Component E 2 of the electric field and component ε 23 of the strain field within the twophase laminate composite. Unitary E 2 and ε 23 imposed.

FIG. 2 :

 2 FIG.2: Component E 1 of the electric field and component ε 13 of the strain field within the fiber composite with square arrangement for an unitary uniaxial applied electric field E 1 .

TABLE I :

 I Symbols and units

	Physical quantity	Symbol	Tensor rank	Units
	Stress

TABLE II :

 II Material properties of the constituents (After Pettermann and Suresh 2 ). Elastic moduli are in GPa, dielectric moduli are in nC/Vm and piezoelectric moduli are in C/m 2 .

		C11	C12	C13	C33 C44 /2 C66 /2	ȇ31	ȇ33	ȇ15 / √	2	γ 11	γ 33
	Epoxy	8.	4.4	4.4	8.	1.8	1.8	0.	0.	0.		3.72e-02 3.72e-02
	PZT 154.837 83.237 82.712 131.39 25.696 35.8 -2.120582 9.52183 9.34959 4.065	2.079
	A. Laminate composite

TABLE III :

 III Comparison of FFT and FEM results 2,3 for the effective moduli of the fibrous composite with 60% volume fraction of fibers. Units are the same as TableII.

			FFT	FEM 2	FEM 3
	C11		25.11	25.16	25.14
	C12		8.57	8.73	8.68
	C13		11.78	11.85	11.83
	C33		54.46	54.61	54.52
	C44 /2		6.63	6.68	6.64
	C66 /2		4.62	4.64	4.64
	ȇ31		-0.20	-0.20	-0.20
	ȇ33		6.43	6.45	6.45
	ȇ15 / √	2	0.052	0.052	0.052
	γ 11		0.157	0.158	0.157
	γ 33		1.277	1.282	1.280

  The coupled constitutive equations can thus be expressed in the following concise formJ α = L αβ F β , (α, β) ∈ [1, 9] 2 , (A2) with [J ] = σ 11 σ 22 σ 33 √ 2σ 23 √ 2σ 13 √ 2σ 12 D 1 D 2 D 3 , [F ] = ε 11 ε 22 ε 33 √ 2ε 23 √ 2ε 13 √ 2ε 12 E 1 E 2 E 3 ,

	A1) where (α, β) ∈ [1, 6] 2 and (i, j) ∈ [1, 3] 2 . In the orthonormal basis (e i , i ∈ [1, 3]), the components of σ, ε, C and ȇ are [ σ] = σ 11 σ 22 σ 33 √ 2σ 23 √ 2σ 13 √ 2σ 12 , [ ε] = ε 11 ε 22 ε 33 √ 2ε 23 √ 2ε 13 √ 2ε 12 , [L] =                       C 1111 C 1122 C 1133 √ 2C 1132 √ 2C 1131 √ 2C 1121 -e 111 -e 211 -e 311 C 2211 C 2222 C 2233 √ 2C 2232 √ 2C 2231 √ 2C 2221 -e 122 -e 222 -e 322 C 3311 C 3322 C 3333 √ 2C 3332 √ 2C 3331 √ 2C 3321 -e 133 -e 233 -e 333 √ 2C 2311 √ 2C 2322 √ 2C 2333 2C 2332 2C 2331 2C 2321 -√ 2e 123 -√ 2e 223 -√ 2e 323 √ 2C 1311 √ 2C 1322 √ 2C 1333 2C 1332 2C 1331 2C 1321 -√ 2e 113 -√ 2e 213 -√ √ 2C 1211 √ 2C 1222 √ 2C 1233 2C 1232 2C 1231 2C 1221 -√ 2e 112 -√ 2e 212 -e 111 e 122 e 133 √ 2e 123 √ 2e 113 2e 112 γ 11 γ 12 e 211 e 222 e 233 √ 2e 223 √ 2e 213 2e 212 γ 12 γ 22 e 311 e 322 e 333 √ 2e 323 √ 2e 313 2e 312 γ 13 γ 23 γ 33 √ γ 23 √ γ 13 √ 2e 312 √ 2e 313	                     	.
	C =	            	C 1111 C 2211 C 3311 √ 2C 2311 √ 2C 1311 √ 2C 1211	C 1122 C 2222 C 3322 √ 2C 2322 √ 2C 1322 √ 2C 1222	C 1133 C 2233 C 3333 √ 2C 1233 2C 1232 2C 1231 2C 1221 √ 2C 1132 √ 2C 1131 √ 2C 1121 √ 2C 2232 √ 2C 2231 √ 2C 2221 √ 2C 3332 √ 2C 3331 √ 2C 3321 √ 2C 1333 2C 1332 2C 1331 2C 1321 √ 2C 2333 2C 2332 2C 2331 2C 2321	            	,
	[ȇ] =	    	e 111 e 122 e 133 e 211 e 222 e 233 e 311 e 322 e 333	√ √ √	2e 123 2e 223 2e 323	√ √ √	2e 113 2e 213 2e 313	√ √ √	2e 112 2e 212 2e 312	    	.

Appendix B: Expression of the influence factors for the binary laminate composite

Analytic expressions of the influence factors of relation ( 12) can be found in 22 . They read, for phase (1),

) .

(B1)


Anne-Claire Egloffe

Lipschitz stability estimate for the Stokes system with Robin boundary conditions

Keywords: Inverse problem, Lipschitz stability estimate, Stokes system Problème inverse, Inégalité de stabilité Lipschitzienne, Système de Stokes

In this research report, we study the inverse problem of identifying a Robin coefficient defined on some non accessible part of the boundary from measurements available on the other part of the boundary, for (u, p) solution of the Stokes system. We prove a Lipschitz stability estimate, under the a priori assumption that the Robin coefficient is piecewise constant. To do so, we use unique continuation estimates for the Stokes system proved in [BEG12] and the approach developed by E. Sincich in [Sin07] to solve a similar inverse problem for the Laplace equation.

Introduction

Let Ω ⊂ R d , with d = 2, 3, be a Lipschitz bounded connected open set such that ∂Ω = Γ l ∪ Γ 0 ∪ Γ out .

We are interested in the following system:

                 -∆u + ∇p = 0, in Ω, div u = 0, in Ω, u = 0, in Γ l , ∂u ∂n -pn = g, on Γ 0 , ∂u ∂n -pn + qu = 0, on Γ out , (1)
where Γ out = N i=1 Γ i and n the exterior normal to Ω. The aim of this paper is to identify the Robin coefficient q defined on the non accessible part of the boundary Γ out from available data on Γ 0 . We provide a Lipschitz stability estimate, under the a priori assumption that the Robin coefficient is piecewise constant on Γ out :

q |Γi = q i ∈ R + , for 1 ≤ i ≤ N.
(2)

Such problems can be viewed as a generalization of some problems which appear naturally in the modeling of biological problems like, for example, blood flow in the cardiovascular system (see [START_REF] Quarteroni | Analysis of a geometrical multiscale model based on the coupling of ODEs and PDEs for blood flow simulations[END_REF] and [START_REF] Vignon-Clementel | Outflow boundary conditions for three-dimensional finite element modeling of blood flow and pressure in arteries[END_REF]) or airflow in the lungs (see [START_REF] Baffico | Multiscale modeling of the respiratory tract[END_REF]). For an introduction on the modeling of the airflow in the lungs and on different boundary conditions which may be prescribed, we refer to [START_REF] Egloffe | Étude de quelques problèmes inverses pour le système de Stokes[END_REF]. The part of the boundary Γ 0 represents a physical boundary on which measurements are available and Γ out represents artificial boundaries on which Robin boundary conditions or mixed boundary conditions involving the fluid stress tensor and its flux at the outlet are prescribed.

The uniqueness for such problems is not an issue and has already been investigate in [START_REF] Boulakia | Stability estimates for a robin coefficient in the two-dimensional stokes system[END_REF]. It is obtained as a corollary of a unique continuation result for the Stokes system proven by C. Fabre and G. Lebeau (see [START_REF] Fabre | Prolongement unique des solutions de l'equation de Stokes[END_REF]). Concerning the stability, logarithmic stability estimates have been obtained in [START_REF] Boulakia | Stability estimates for a robin coefficient in the two-dimensional stokes system[END_REF] and [START_REF] Boulakia | Unique continuation estimates for the stokes system; application to an inverse problem[END_REF] for more general than piecewise constant Robin coefficients. The main tools used in both cases are Carleman inequalities, global in the first case and local in the second one. Let us point out that, due to the mixed boundary conditions, the solution of system (1) is not regular in a neighborhood of the junction between two different boundary conditions. Thus, we can not expect in general that the solution belongs, at least, to H 2 (Ω) × H 1 (Ω). As a consequence, we can not use global Carleman inequalities requiring global regularity on the solution. The Lipschitz stability estimate obtained in this paper is based on local regularity on the solution of system (1) and on the open set Ω and on the stability estimates for the unique continuation property of the Stokes system proved in [START_REF] Boulakia | Unique continuation estimates for the stokes system; application to an inverse problem[END_REF].

Comparable problems have been widely studied for the Laplace equation (see for instance [START_REF] Alessandrini | Stable determination of corrosion by a single electrostatic boundary measurement[END_REF], [START_REF] Bellassoued | Stability estimate for an inverse boundary coefficient problem in thermal imaging[END_REF], [START_REF] Chaabane | Logarithmic stability estimates for a Robin coefficient in two-dimensional Laplace inverse problems[END_REF], [START_REF] Chaabane | Identification of Robin coefficients by the means of boundary measurements[END_REF], [START_REF] Cheng | Stable determination of a boundary coefficient in an elliptic equation[END_REF] and [START_REF] Sincich | Lipschitz stability for the inverse Robin problem[END_REF]. In this case, it is in general a problem arising in corrosion detection which consists of determining a Robin coefficient on the inaccessible portion of the boundary by electrostatic measurements performed on the accessible one. Most of the paper obtained logarithmic stability estimates but under some restricting assumption on the Robin coefficient and on the flux g, it is possible to obtain Lipschitz stability estimate. For instance, S. Chaabane and M. Jaoua obtained in [START_REF] Chaabane | Identification of Robin coefficients by the means of boundary measurements[END_REF] both local and monotone global Lipschitz stability for regular Robin coefficient and under the assumption that the flux g is non negative. Relaxing this constraint, they obtained in [START_REF] Alessandrini | Stable determination of corrosion by a single electrostatic boundary measurement[END_REF] a logarithmic stability estimate. More recently, E. Sincich has obtained in [START_REF] Sincich | Lipschitz stability for the inverse Robin problem[END_REF] a Lipschitz stability estimate under the further a priori assumption of a piecewise constant Robin coefficient.

Let us explain the structure of this paper. In Section 2, we begin with existence and regularity results on (u, p) solution of system (1): we state global and local regularity results which will be useful to prove the Lipschitz stability estimate. Then, we give in Section 3 technical lemmas and unique continuation estimates which will be useful to prove the Lipschitz stability estimate. Finally, Section 4 is dedicated to the main result: we state and prove the Lipschitz stability estimate.

In the following, we will not distinguish vector valued functions and scalar valued functions. Moreover, when we are not more specific, C > 0 is a constant whose value may change from a line to an other. Let us introduce some notations that we will be useful throughout this paper.

Notation 1.1. For x ∈ R d and r > 0, we denote B r (x) the ball of center x and of radius r.

Notation 1.2. Let Γ ⊂ ∂Ω be a non empty part of the boundary. We denote by

Γ in = {x ∈ Γ/d(x, ∂Ω\Γ) > 0}.
Since the open set Ω is Lipschitz, it satisfies the cone property: Definition 1.3 (Cone property). We say that Ω satisfies the cone property if there exists θ ∈ (0, π 2 ) and R 0 > 0 such that for all x 0 ∈ ∂Ω, there exists ξ ∈ R d , |ξ| = 1 such that the finite cone

C = {x ∈ R d /(x -x 0 ) • ξ > |x -x 0 | cos θ, |x -x 0 | ≤ R 0 } is included in Ω.

Regularity results

We focus in this section on global and local regularity results for system (1). Even if global H 2 ×H 1 regularity is not expected in general due to the mixed Dirichlet and Neumann boundary conditions, we obtain local regularity results inside the domain and near the boundary, as long as we stay away from the junction between two different boundary conditions. We need to introduce functional spaces:

V Γ l = {v ∈ H 1 (Ω)/v |Γ l = 0 and div v = 0}, (3)
and

H Γ l = V L 2 (Ω) Γ l . (4) Moreover, for g ∈ H -1 2 (∂Ω) and v ∈ H 1 2 (∂Ω), we denote by < g, v > -1 2 , 1
2 ,∂Ω the image of v by the linear form g.

Global regularity

Proposition 2.1. Let R M > 0, 1 Γ0 g ∈ H -1 2 (∂Ω) and assume that q satisfies (2). Then, system (1) admits a unique solution (u, p) ∈ V Γ l × L 2 (Ω). Moreover, if we assume the q ≤ R M , there exists a constant C(R M ) > 0 such that:

u H 1 (Ω) + p L 2 (Ω) ≤ C(R M ) 1 Γ0 g H -1 2 (∂Ω)
.

(5)

Inria

Proof of Proposition 2.1. The variational formulation of the problem is: find u ∈ V Γ l such that for every v ∈ V Γ l ,

Ω ∇u : ∇v + Γout qu • v =< 1 Γ0 g, v > -1 2 , 1 2 ,∂Ω . (6)
We denote by

a q (u, v) = Ω ∇u : ∇v + Γout qu • v, (7)
and

L(v) =< 1 Γ0 g, v > -1 2 , 1
2 ,∂Ω . We easily verify that a q is a continuous symmetric bilinear form on V Γ l . According to Poincaré inequality, the bilinear form a q is coercive on V Γ l . On the other hand, L is a continuous linear form on V Γ l . Thus we get the existence and uniqueness of u ∈ V Γ l solution of equations (1) using the Lax-Milgram Theorem.We prove the existence and uniqueness of p ∈ L 2 (Ω) in a classical way, by using De Rham Theorem and the Neumann boundary conditions.

Local regularity

Inside the domain Ω, we have local regularity: this is resumed in Proposition 2.2. Moreover, locally near the boundary, as long as we stay away from the junction between two different boundary conditions, we can also obtain local regularity. We refer to Propositions 2.3 and 2.5 for a statement of these regularity results.

Proposition 2.2. Let R M > 0, m ∈ N * and ω ⊂ Ω be a relatively compact open set. Let 1 Γ0 g ∈ H -1
2 (∂Ω) and assume that q satisfies (2). Then, the solution (u, p) of system (1) belongs to H m+1 (ω)×H m (ω). Moreover, if we assume that q ≤ R M , there exists a constant C(R M ) > 0 such that:

u H m+1 (ω) + p H m (ω) ≤ C(R M ) 1 Γ0 g H -1 2 (∂Ω) . ( 8
)
Although the proof of this result is classical, we give here a sketch of the proof for the sake of completeness.

Proof of Proposition 2.2. We prove this proposition by induction on m. For m = 0, the result is given by Proposition 2.1.

We assume that Proposition 2.2 holds for some fixed m. Let us prove that the proposition also holds for m + 1. Let ω be an open set of class C m+1,1 such that ω ⊂ ω ⊂ Ω. We localize in the neighborhood of ω. Let χ ∈ C ∞ c (ω) be such that χ = 1 in ω and 0 ≤ χ ≤ 1 everywhere else. We denote by (v, π) = (χu, χp). Note that (v, π) is solution if the following problem:

   -∆v + ∇π = -∆χu -2∇u∇χ + ∇χp, in ω, div v = ∇χ • u, in ω, ∂v ∂n -πn = 0, on ∂ω.
Let us denote by f = -∆χu -2∇u∇χ + ∇χp and h = ∇χ • u. By the induction assumption, we deduce that the (f, h) belongs to H m (ω) × H m+1 (ω). By application of regularity result for the Stokes system with Neumann boundary condition (see [START_REF] Boyer | Éléments d'analyse pour l'étude de quelques modèles d'écoulements de fluides visqueux incompressibles[END_REF]), we deduce that (v, π) belongs to H m+2 (ω) × H m+1 (ω). Since χ = 1 on ω, we obtain the desired result.

We now study the regularity near the boundary of the domain. Proposition 2.3 states regularity result in the restriction to Ω of a neighborhood of any point x j ∈ Γ in j , for i = j, . . . , N .

Proposition 2.3. Let m ∈ N * , R M > 0, 1 ≤ j ≤ N and x j ∈ Γ in j . We assume that Γ j is of class C m,1 if m ≥ 1 and Lipschitz otherwise. Let 1 Γ0 g ∈ H -1
2 (∂Ω) and assume that q satisfies (2). Then, there exists R > 0 such that the solution (u, p) of system (1) belongs to

H m+1 (B R (x j )∩ Ω)×H m (B R (x j )∩Ω). Furthermore, if we assume that q ≤ R M , there exists a constant C(R M ) > 0 such that u H m+1 (B R (xj )∩Ω) + p H m (B R (xj )∩Ω) ≤ C(R M ) 1 Γ0 g H -1 2 (∂Ω)
.

The proof of Proposition 2.3 relies on classical argument and is similar to the proof of Proposition 2.2.

Then, we deduce from Proposition 2.3 local Hölder regularity on (u, p) solution of system (1) near the boundary Γ j , as long as we stay away from the junction between two different boundary conditions.

Corollary 2.4. Let 1 ≤ j ≤ N and x j ∈ Γ in j . We assume that Γ j is of class C 2,1 . Let 1 Γ0 g ∈ H -1
2 (∂Ω) and assume that q satisfies (2). Then, there exists R > 0 and 0 < β < 1 such that the solution (u, p) of system (1) belongs to

C 1,β (B R (x j ) ∩ Ω) × C 0,β (B R (x j ) ∩ Ω) for all j = 1, . . . , N . Proof of Corollary 2.4. Let 1 ≤ j ≤ N . Thanks to Proposition 2.3, there exists R j > 0 such that (u, p) ∈ H 3 (B Rj (x j ) ∩ Ω) × H 2 (B Rj (x j ) ∩ Ω). Thanks to Proposition 2.2, we know that there exists a connected open set Ω ⊂ Ω of class C 2,1 such that N j=1 B Rj (x j ) ∩ Ω ⊂ Ω and (u, p) ∈ H 3 ( Ω) × H 2 ( Ω). Then, since for all m ∈ N such that 2(m -1) ≤ d < 2m, there exists 0 < λ < 1 such that H m ( Ω) → C 0,λ ( Ω) (see [Ada75]), we deduce that there exits 0 < β < 1 such that H 2 ( Ω) ⊂ C 0,β ( Ω), which implies that (u, p) ∈ C 1,β (B R (x j ) ∩ Ω) × C 0,β (B R (x j ) ∩ Ω) for all j = 1, . . . , N , with R = min 1≤j≤N R j .
Proposition 2.5 below states regularity result in the restriction to Ω of a neighborhood of any point

x 0 ∈ Γ in 0 . Proposition 2.5. Let R M > 0, m ∈ N and x 0 ∈ Γ in 0 . We assume that Γ 0 is of class C m,1 if m ≥ 1 and Lipschitz otherwise. Let 1 Γ0 g ∈ H m-1
2 (∂Ω) and assume that q satisfies (2). Then, there exists R > 0 such that the solution (u, p) of system (1) belongs to

H m+1 (B R (x 0 )∩ Ω)×H m (B R (x 0 )∩Ω). Furthermore, if we assume that q ≤ R M , there exists a constant C(R M ) > 0 such that u H m+1 (B R (x0)∩Ω) + p H m (B R (x0)∩Ω) ≤ C(R M ) 1 Γ0 g H m-1 2 (∂Ω)
.

3 Preliminary results

Useful lemmas

The following lemmas will be useful throughout this paper.

Lemma 3.1. Let A > 0, B > 0, C 1 > 0, C 2 > 0 and D > 0. We assume that there exists

γ 1 > 0 such that D ≤ Ae C1γ + Be -C2γ , (9)
for all γ ≥ γ 1 and c 0 > 0 such that D ≤ c 0 B. Then, there exists C > 0 such that:

D ≤ CA C 2 C 1 +C 2 B C 1 C 1 +C 2 .

Inria

We refer to [START_REF] Robbiano | Fonction de coût et contrôle des solutions des équations hyperboliques[END_REF] for a proof of this lemma.

Lemma 3.2. Let A ∈ R, µ ∈ R * and (β k ) k∈N ∈ R N .
If, for all k ∈ N * , we have

β k ≤ 1 µ k-1 (β k-1 ) α A 1-α , ( 10
)
then β k ≤ 1 µ ι β α k 0 A 1-α k ,
where ι = k-1 j=1 jα k-1-j . This lemma is proved in [START_REF] Bourgeois | About stability and regularization of ill-posed elliptic Cauchy problems: the case of Lipschitz domains[END_REF]. For the sake of completeness, we write it.

Proof of Lemma 3.2. We rewrite inequality (10) under the form:

β k A ≤ 1 µ k-1 β k-1 A α .
By iterating the above inequality, we get:

β k A ≤ 1 µ ι β 0 A α k , where ι = k-1 j=1 jα k-1-j .
Lemma 3.3. Let A, B, C 1 and D be positive numbers and 0 < α < 1. Assume that

D ≤ C 1 A α B 1-α .
Then, for all > 0 D ≤ c A + s B,

where c = C 1 α 1 and s = α 1 -α . Proof of Lemma 3.3. Let > 0. We rewrite C 1 A α B 1-α = C 1/α 1 A α α B 1-α .
Then it is sufficient to apply Young inequality:

D ≤ α C 1/α 1 A + (1 -α) s B.
Since 0 < α < 1, the desired inequality follows.

Unique continuation estimates

In this section, we state some unique continuation estimates for the Stokes system which will be useful in the next subsection to prove the Lipschitz stability estimate. They are obtained as corollaries of unique continuation estimates proved in [START_REF] Boulakia | Unique continuation estimates for the stokes system; application to an inverse problem[END_REF]. Let us begin by recalling a proposition from [START_REF] Boulakia | Unique continuation estimates for the stokes system; application to an inverse problem[END_REF] which allows to transmit information from a part of the boundary of Ω to a relatively compact open set in Ω.

Proposition 3.4. Assume that D is of class C ∞ . Let 0 < ν ≤ 1 2 .
Let Γ be a non empty open subset of the boundary of D. Let ω be a relatively compact open set in D. Then, there exists C, σ > 0, such that for all > 0 and for all (u, p) ∈ H

3 2 +ν (D) × H 3 2 +ν (D) solution of -∆u + ∇p = 0, in D, div u = 0, in D, (11)
u H 1 (ω) + p H 1 (ω) ≤ C u H 1 (Γ) + p H 1 (Γ) + ∂u ∂n L 2 (Γ) + ∂p ∂n L 2 (Γ) + σ ( u H 1 (D) + p H 1 (D) ).
In our case, we can not apply directly Proposition 3.4 in the open set Ω to (u, p) solution of system (1) because regularity is needed both on the solution of the Stokes system (11) and on the open set. Nevertheless, Proposition 3.5 below is obtain as a corollary of Proposition 3.4 and allows us to transmit information from a part of the boundary Γ ⊂ Γ 0 to a relatively compact open ω set included in Ω.

Proposition 3.5. Assume that Γ 0 is of class C ∞ . Let R M > 0, M 1 > 0, Γ ⊆ Γ 0 be a non empty open subset of the boundary of Ω such that Γ ∩ Γ l Γ ∩ Γ out = ∅ and ω ⊂ Ω be a relatively compact open set. Let 1 Γ0 g ∈ H 3 2 (∂Ω) be such that 1 Γ0 g H 3 2 (∂Ω)
≤ M 1 and assume that q satisfies (2) and q ≤ R M .

Then, there exist constants C(R M , M 1 ) > 0 and 0 < δ < 1 such that for all (u, p) solution of system (1), the following inequality is satisfied:

u H 1 (ω) + p H 1 (ω) ≤ C(R M , M 1 ) u L 2 (Γ) + p L 2 (Γ) + ∂u ∂n L 2 (Γ) + ∂p ∂n L 2 (Γ) δ . (12)
Proof of Proposition 3.5. Thanks to local regularity results stated in Subsection 2.2, we know that there exists a connected open set Ω ⊂ Ω of class C ∞ such that Γ ⊂ ∂ Ω, ω ⊂ Ω and such that the solution (u, p) of (1) belongs to H 3 ( Ω) × H 2 ( Ω). Moreover, there exists a constant C(R M , M 1 ) > 0 such that:

u H 3 ( Ω) + p H 2 ( Ω) ≤ C(R M , M 1 ). (13)
We apply Proposition 3.4: there exists σ > 0 and C > 0 such that, for all ˜ > 0,

u H 1 (ω) + p H 1 (ω) ≤ C ˜ u H 1 (Γ) + p H 1 (Γ) + ∂u ∂n L 2 (Γ) + ∂p ∂n L 2 (Γ) + ˜ σ u H 1 ( Ω) + p H 1 ( Ω) . ( 14
)
Note that it is the H 1 norms of u and p on Γ which appear in the first term in the right hand-side of (14). In order to replace them with the L 2 norms of u and p on Γ, we use an interpolation inequality: there exists c > 0 such that

u H 1 (Γ) + p H 1 (Γ) ≤ c u 1 3 L 2 (Γ) u 2 3 H 3 2 (Γ) + p 1 3 L 2 (Γ) p 2 3 H 3 2 (Γ)
.

Let > 0. If we write

u 1 3 L 2 (Γ) u 2 3 H 3 2 (Γ) = 1 u L 2 (Γ) 1 3 1 2 u H 3 2 (Γ) 2 3 , Inria and p 1 3 L 2 (Γ) p 2 3 H 3 2 (Γ) = 1 p L 2 (Γ) 1 3 1 2 p H 3 2 (Γ) 2 3 ,
according to Young inequality and to the trace injection H 2 ( Ω) → H 3 2 (Γ), we obtain:

u H 1 (Γ) + p H 1 (Γ) ≤ c 1 2 u H 3 ( Ω) + p H 2 ( Ω) + 1 u L 2 (Γ) + p L 2 (Γ) . (15)
Let > 0. By combining inequalities (15) with = ˜ 2(σ+1) and inequalities ( 14) with ˜ = 1 2σ+3 , we obtain the existence of C > 0 and s > 0 such that for all > 0:

u H 1 (ω) + p H 1 (ω) ≤ C 1 u L 2 (Γ) + p L 2 (Γ) + ∂u ∂n L 2 (Γ) + ∂p ∂n L 2 (Γ) + s u H 3 ( Ω) + p H 2 ( Ω) . (16)
We conclude by using Lemma 3.1 and inequality (13).

Then, Lemma 3.6 allows to transmit information inside the domain Ω, from a relatively compact open set include in Ω to an other. Lemma 3.6. [Three balls inequality] Let ρ > 0 and q ∈ R d . There exist C > 0, 0 < α < 1 such that for all function (u, p) ∈ H 1 (B 8ρ (q)) × H 1 (B 8ρ (q)) solution of

-∆u + ∇p = 0, div u = 0, (17)
in B 8ρ (q) the following inequality is satisfied:

u H 1 (B3ρ(q)) + p L 2 (B3ρ(q)) ≤ C u H 1 (Bρ(q)) + p L 2 (Bρ(q)) α u H 1 (B8ρ(q)) + p L 2 (B8ρ(q)) 1-α , ( 18
)
with α = g(3ρ) -g 7 2 ρ g( ρ 3 ) -g 7 2 ρ
and with g(r) = e -λr 2 , for λ large enough.

In the following, in order to refer to this inequality, we will say that the three balls inequality associated to q and ρ is satisfied, with the associated constants C > 0 and α > 0.

Remark 3.7. Lemma 3.6 is the counterpart, in the case of the Stokes system, of the so-called three balls inequality for the Laplacian. We refer to [START_REF] Bourgeois | About stability and regularization of ill-posed elliptic Cauchy problems: the case of Lipschitz domains[END_REF] or [START_REF] Lin | Quantitative uniqueness for the power of laplacian with singular coefficients[END_REF] for a three balls inequality for the Laplacian. Note that in [START_REF] Lin | Optimal three-ball inequalities and quantitative uniqueness for the Stokes system[END_REF], C.-H. Lin, G. Uhlmann and J.-N. Wang have obtained an optimal three balls inequality for the Stokes system involving only the velocity in the L 2 norm. From this inequality, they derive an upper bound on the vanishing order of any non trivial solution u to the Stokes system.

The ideas of the proof are the same as those developed in [START_REF] Boulakia | Unique continuation estimates for the stokes system; application to an inverse problem[END_REF]. We refer to [START_REF] Egloffe | Étude de quelques problèmes inverses pour le système de Stokes[END_REF] for a complete proof of this result. In the proof of the Lipschitz stability estimate, we will apply Lemma 3.6 to a sequence of balls with decreasing radius which approachs the boundary. To this aim, we need to know the behavior of the constants when we pass from a ball to another. This is done in Lemma 3.8 below.

Lemma 3.8. Let ρ > 0, (q, q) ∈ R d × R d , and µ ∈ (0, 1). We denote by ρ = µρ. We assume that the three balls inequality (18) associated to q and ρ holds for some constants C > 0 and α > 0.

Then, for all functions (u, p) ∈ H 1 (B 8ρ (q)) × H 1 (B 8ρ (q)) solution of (17) in B 8ρ (q), the following inequality is satisfied:

u H 1 (B3ρ(q)) + p L 2 (B3ρ(q)) ≤ C u H 1 (Bρ(q)) + p L 2 (Bρ(q)) α u H 1 (B8ρ(q)) + p L 2 (B8ρ(q)) 1-α , ( 19
)
where C = C µ > 0. In other words, the three balls inequality associated to q and ρ is satisfied with the associated constants C µ > 0 and α > 0.

Proof of Lemma 3.8. This lemma is inspired from [START_REF] Bourgeois | About stability and regularization of ill-posed elliptic Cauchy problems: the case of Lipschitz domains[END_REF] where L. Bourgeois and J. Dardé are concerned with the operator P k = -∆ -k, with k ∈ R, and use similar techniques. Let m ∈ N * . By performing the change of variables B mρ (q) → B mρ (q), x → q + (x -q)µ, we get:

Bmρ(q) |u(x)| 2 + |∇u(x)| 2 dx = µ d Bmρ(q) |u(q + µ(x -q))| 2 + |∇u(q + µ(x -q))| 2 dx.
Let us denote by u(x) = u(q + µ(x -q)) and p(x) = µp(q + µ(x -q)). Noticing that ∇u(x) = µ∇u(q + µ(x -q)), we get:

Bmρ(q) |u(x)| 2 + |∇u(x)| 2 dx = µ d Bmρ(q) |u(x)| 2 dx + 1 µ 2 |∇u(x)| 2 dx.

Moreover, since

Bmρ(q) |p(x)| 2 dx = µ d-2 Bmρ(q) |p(x)| 2 dx,
and since 0 < µ < 1, we obtain:

µ d 2 u H 1 (Bmρ(q)) + p L 2 (Bmρ(q)) ≤ u H 1 (Bmρ(q)) + p L 2 (Bmρ(q)) ≤ µ d 2 -1 u H 1 (Bmρ(q)) + p L 2 (Bmρ(q)) .
Observe that (u, p) is solution in B 8ρ (q) of system (17):

-∆u(x) + ∇p(x) = µ 2 (-∆u(x) + ∇p(x)) = 0, for x ∈ B 8ρ (q) and where x = q + (x -q)µ ∈ B 8ρ (q). Thus, (u, p) satisfies (18) for C > 0 and α > 0. We deduce that:

u H 1 (B3ρ(q)) + p L 2 (B3ρ(q)) ≤ µ d 2 -1 u H 1 (B3ρ(q)) + p L 2 (B3ρ(q)) ≤ µ d 2 -1 C u H 1 (Bρ(q)) + p L 2 (Bρ(q)) α u H 1 (B8ρ(q)) + p L 2 (B8ρ(q)) 1-α ≤ C µ u H 1 (Bρ(q)) + p L 2 (Bρ(q)) α u H 1 (B8ρ(q)) + p L 2 (B8ρ(q)) 1-α . Inria 4 Main result Theorem 4.1. Assume that Γ 0 is of class C ∞ and Γ i is of class C 2,1 for i = 1, . . . , N . Let m > 0, R M > 0, M 1 > 0, Γ ⊆ Γ 0 be a non empty open subset of the boundary of Ω such that Γ ∩ Γ l Γ ∩ Γ out = ∅ and let g ∈ H 3
2 (Γ 0 ) be non identically zero on Γ 0 and such that g

H 3 2 (Γ0)
≤ M 1 . We assume that q k satisfies (2) with q i = q k i be such that q k i ≤ R M for i = 1, . . . , N and k = 1, 2. Let us denote by (u k , p k ) the solution of (1) associated to q = q k for k=1,2. We assume that there exists x j ∈ Γ in j such that |u 2 (x j )| > m, for all j = 1, . . . , N . Then, there exists C(R M , M 1 , N, m) > 0 such that

q 1 -q 2 L ∞ (Γout) ≤ C(R M , M 1 , N, m) u 1 -u 2 L 2 (Γ) + p 1 -p 2 L 2 (Γ) + ∂p 1 ∂n - ∂p 2 ∂n L 2 (Γ) . ( 20
)
Remark 4.2. Since g is not identically zero on Γ 0 , we know, thanks to the uniqueness result (see [START_REF] Boulakia | Stability estimates for a robin coefficient in the two-dimensional stokes system[END_REF]), that for all j = 1, . . . , N , there exists x j ∈ Γ in j such that u 2 (x j ) = 0. We notice however that the constant involved in the estimate (20) depends on u 2 through the constant m. Finding a uniform lower bound for a solution u of system (1) remains an open question.

Remark 4.3. Let 1 ≤ i ≤ N and x i ∈ Γ in i be such that |u 2 (x i )| > m.
Let us give an idea of how information goes from Γ to a neighborhood V i of x i . Thanks to Proposition 3.5, information goes from Γ to a relatively compact open set B 0 included in Ω. Then, we use a sequence of balls (B k ) k∈N with decreasing radius to approach the boundary near x i , taking into account Lemma 3.8. Finally, we use the boundary condition on Γ i and the local Hölder regularity of the solution on V i (see Corollary 2.4). We refer to Figure 1 for an illustration.

Proof of Theorem 4.1. We follow the approach developed in [START_REF] Sincich | Lipschitz stability for the inverse Robin problem[END_REF] in the case of the Laplace equation.

We consider:

(w, π) = u 1 -u 2 N j=1 |q 1 j -q 2 j | , p 1 -p 2 N j=1 |q 1 j -q 2 j | . ( 21
)
According to Proposition 2.1, (w, π) belongs to V Γ l × L 2 (Ω). Since for k = 1, 2, q k is piecewise constant, (w, π) is solution of:

                   -∆w + ∇π = 0, in Ω, div w = 0, in Ω, w = 0, on Γ l , ∂w ∂n -πn = 0, on Γ 0 , ∂w ∂n -πn + q 1 w = (q 2 -q 1 ) N j=1 |q 1 j -q 2 j | u 2 , on Γ out . ( 22
)
Observe that if we do not assume that q k is piecewise constant, we get additional terms which depend on the derivative of q k in system (22).

Step 1 : Since the open set Ω satisfies the cone property, there exists θ ∈ (0, 1) and R 0 > 0 such that for all i = 1, . . . , N there exists

ξ i ∈ R d , |ξ i | = 1 such that the finite cone C i = {x ∈ R d /(x -x i ) • ξ i > |x -x i | cos θ and |x -x i | ≤ R 0 } is included in Ω.
We are going to construct a sequence of balls (B ρ k (ζ i,k )) k∈N with decreasing radius and whose center is converging through x i . For this sequence of balls, we will prove that there exists 0 < α < 1, 0 < µ < 1 and C > 0 such that for all i = 1, . . . , N and k ∈ N, the following estimate is satisfied for all > 0:

Γ x i B 0 B 1 V x i
w H 1 (Bρ k (ζ i,k )) + π L 2 (Bρ k (ζ i,k )) ≤ e C α k log 1 µ k-1 w H 1 (Bρ 0 (ζi,0)) + π L 2 (Bρ 0 (ζi,0)) + C α k w H 3 ( Ωi) + π H 2 ( Ωi) , (23)
where Ωi ⊂ Ω is an open set such that Γ ⊂ ∂ Ωi and C i ⊂ Ωi and such that (w, π) ∈ H 3 ( Ωi ) × H 2 ( Ωi ). We consider

C i = {x ∈ R d /(x -x i ) • ξ i > |x -x i | cos θ et |x -x i | ≤ R 0 }, with θ = arcsin(t sin θ). (24)
The parameter t belongs to (0, 1) and will be specified later on. Note that we have

C i ⊂ C i ⊂ Ω.
We denote by

ζ i,0 = x i + R 0 2 ξ i , d 0 = |ζ i,0 -x i | = R 0 2
and ρ 0 = d 0 sin θ . For k ∈ N * , we define the sequence of balls by induction:

ζ i,k+1 = x i + µ|ζ i,k -x i |ξ i , d k+1 = |ζ i,k+1 -x i |, ρ k+1 = d k+1 sin(θ ), with µ = 1 -sin θ 1 + sin θ ⇐⇒ sin θ = 1 -µ 1 + µ . (25)
We refer to Figure 2 for an illustration of this construction. This construction implies that d k+1 = µd k and ρ k+1 = µρ k . We choose the parameter t involved in (24) such that Next, we apply Lemma 3.6 with ρ = ρ 0 and q = ζ i,0 : there exists C > 0 and α > 0 such that

B 8ρ k (ζ i,k ) ⊂ C i ⊂ Ω for all k ∈ N, that is to say, using (24), 8ρ k = 8d k sin θ ≤ sin θd k ⇔ t ≤ 1 8 . Inria ζ i,k x i ϑ' ϑ d k ρ k ξ i ! ✗ ζ i,0
w H 1 (B3ρ 0 (ζi,0)) + π L 2 (B3ρ 0 (ζi,0)) ≤ C w H 1 (Bρ 0 (ζi,0)) + π L 2 (Bρ 0 (ζi,0)) α w H 1 (B8ρ 0 (ζi,0)) + π L 2 (B8ρ 0 (ζi,0)) 1-α .
Taking into account Lemma 3.8, since ρ k-1 = µ k-1 ρ 0 , it follows that, for all k ∈ N * :

w H 1 (B 3ρ k-1 (ζ i,k-1 )) + π L 2 (B 3ρ k-1 (ζ i,k-1 )) ≤ C µ k-1 w H 1 (Bρ k-1 (ζ i,k-1 )) + π L 2 (Bρ k-1 (ζ i,k-1 )) α × w H 1 (B 8ρ k-1 (ζ i,k-1 )) + π L 2 (B 8ρ k-1 (ζ i,k-1 )) 1-α . (26)
By construction, we have:

B ρ k (ζ i,k ) ⊆ B 3ρ k-1 (ζ i,k-1 ). (27)
Indeed, we have

ζ i,k -ζ i,k-1 = µ(d k-1 -d k-2 )ξ i = µ sin θ (ρ k-1 -ρ k-2 )ξ i = - µ sin θ (1 -µ)ρ k-2 ξ i and using (24), we obtain ζ i,k -ζ i,k-1 = -µ(1 + µ)ρ k-2 ξ i . Then, if x ∈ B ρ k (ζ i,k ), we have: |x -ζ i,k-1 | ≤ |x -ζ i,k | + |ζ i,k-1 -ζ i,k | ≤ ρ k + µ(1 + µ)ρ k-2 ≤ 3ρ k-1 .
By combining (26) and ( 27) we deduce that there exists C > 0 and 0 < α < 1 such that for all k ∈ N * , we have:

w H 1 (Bρ k (ζ i,k )) + π L 2 (Bρ k (ζ i,k )) ≤ C µ k-1 w H 1 (Bρ k-1 (ζ i,k-1 )) + π L 2 (Bρ k-1 (ζ i,k-1 )) α w H 1 (B 8ρ k-1 (ζ i,k-1 )) + π L 2 (B 8ρ k-1 (ζ i,k-1 )) 1-α . ( 28
)

RR n°8222

Let Ωi ⊂ Ω be an open set such that Γ ⊂ ∂ Ωi , C i ⊂ Ωi and such that (w, π) ∈ H 3 ( Ωi ) × H 2 ( Ωi ). Note that such an open set exists thanks to local regularity results stated in Subsection 2.2. We have, for all k ∈ N * :

w H 1 (B8ρ k-1 (ζ i,k-1 )) + π L 2 (B8ρ k-1 (ζ i,k-1 )) ≤ w H 3 ( Ωi) + π H 2 ( Ωi) .
Thus, we can rewrite inequality (28) as:

w H 1 (Bρ k (ζ i,k )) + π L 2 (Bρ k (ζ i,k )) ≤ C µ k-1 w H 1 (Bρ k-1 (ζ i,k-1 )) + π L 2 (Bρ k-1 (ζ i,k-1 )) α w H 3 ( Ωi) + π H 2 ( Ωi) 1-α .
We now apply Lemma 3.2 with

β k = w H 1 (Bρ k (ζ i,k )) + π L 2 (Bρ k (ζ i,k )) and A = C w H 3 ( Ωi) + π H 2 ( Ωi)
, to obtain:

w H 1 (Bρ k (ζ i,k )) + π L 2 (Bρ k (ζ i,k )) ≤ 1 µ ι w H 1 (Bρ 0 (ζi,0)) + π L 2 (Bρ 0 (ζi,0)) α k C w H 3 ( Ωi) + π H 2 ( Ωi) 1-α k , where ι = k-1 j=1 jα k-1-j . Note that ι ≤ (k -1)   k-2 j=0 α j   ≤ k -1 1 -α .
To summarize, the following inequality is satisfied:

w H 1 (Bρ k (ζ i,k )) + π L 2 (Bρ k (ζ i,k )) ≤ 1 µ k-1 1-α
w H 1 (Bρ 0 (ζi,0)) + π L 2 (Bρ 0 (ζi,0))

α k C( w H 3 ( Ωi) + π H 2 ( Ωi) ) 1-α k . Let > 0. By Lemma 3.3 with A = w H 1 (Bρ 0 (ζi,0)) + π L 2 (Bρ 0 (ζi,0)) , B = C( w H 3 ( Ωi) + π H 2 ( Ωi) ), C 1 = 1 µ k-1 1-α
, and

D = w H 1 (Bρ k (ζ i,k )) + π L 2 (Bρ k (ζ i,k
)) , we obtain:

w H 1 (Bρ k (ζ i,k )) + π L 2 (Bρ k (ζ i,k )) ≤ c k w H 1 (Bρ 0 (ζi,0)) + π L 2 (Bρ 0 (ζi,0)) + C s k w H 3 ( Ωi) + π H 2 ( Ωi) , with s k = α k 1 -α k and c k = 1 µ k-1 1-α 1/α k . Note that s k ≥ α k implies, for 0 < < 1 that: w H 1 (Bρ k (ζ i,k )) + π L 2 (Bρ k (ζ i,k )) ≤ e C α k log 1 µ k-1 w H 1 (Bρ 0 (ζi,0)) + π L 2 (Bρ 0 (ζi,0)) + C α k w H 3 ( Ωi) + π H 2 ( Ωi) .
Let us remark that the previous inequality is clearly satisfied for ≥ 1 by continuity of the injection

H 1 ( Ωi ) → H 1 (B ρ k (ζ i,k ))
. Thus, we obtain inequality (23).

Step 2 : Combining local Hölder regularity of the solution and inequality (23), and then optimizing the resulting inequality, we will obtain the desired Lipschitz stability estimate (20).

Inria

Let us recall that thanks to Corollary 2.4, there exists R > 0 and 0

< β < 1 such that (w, π) ∈ C 1,β (B R (x i ) ∩ Ω) × C 0,β (B R (x i ) ∩ Ω) for all i = 1, . . . , N . Let 0 < < min(R, 1) = 0 . If d k + ρ k < , that is to say µ k (d 0 + ρ 0 ) < , we have B ρ k (ζ i,k ) ⊂ B (x i ).
Let k 0 = k 0 ( ) be the smallest k which satisfies this inequality. We have:

log((d 0 + ρ 0 )/ ) log (1/µ) ≤ k 0 < log((d 0 + ρ 0 )/ ) log (1/µ) + 1. (29)
Since (w, π) is solution of system (22), we obtain, using the boundary condition on Γ i :

|q 1 i -q 2 i | N j=1 |q 1 j -q 2 j | |u 2 (x i )| ≤ ∂w ∂n (x i ) +|π(x i )|+|q 1 i ||w(x i )| ≤ |∇w(x i )|+|π(x i )|+R M |w(x i )|. (30) Let y i ∈ B ρ k 0 (ζ i,k0
). Using the Hölder regularity of w, ∇w and π on B ρ k 0 (ζ i,k0 ) ⊂ B (x i ) ⊂ B R (x i ), we have:

|q 1 i -q 2 i | N j=1 |q 1 j -q 2 j | |u 2 (x i )| ≤ |∇w(y i )| + |π(y i )| + R M |w(y i )| + C|x i -y i | β .
Moreover, since |u 2 (x i )| > m, we obtain that, for all y i in B ρ k 0 (ζ i,k0 ),

|q 1 i -q 2 i | N j=1 |q 1 j -q 2 j | m ≤ |∇w(y i )| + |π(y i )| + R M |w(y i )| + C β .
Let us denote by ω d the volume of the unit ball in R d . By integrating in L 2 norm the previous inequality in B ρ k 0 (ζ i,k0 ), we obtain:

|q 1 i -q 2 i | N j=1 |q 1 j -q 2 j | ≤ C(m, R M ) ω 1 2 d ρ d 2 k0 w H 1 Bρ k 0 (ζi,k 0 ) + π L 2 Bρ k 0 (ζi,k 0 ) + C β .
The previous inequality together with (23) yields to the existence of 0 < α < 1, 0 < µ < 1 such that for all i = 1, . . . , N the following estimate holds for all > 0

|q 1 i -q 2 i | N j=1 |q 1 j -q 2 j | ≤ C(m, R M ) ρ d 2 k0 e C α k 0 log 1 µ k 0 -1 w H 1 (Bρ 0 (ζi,0)) + π L 2 (Bρ 0 (ζi,0)) + C(m, R M ) ρ d 2 k0 α k 0 w H 3 ( Ωi) + π H 2 ( Ωi) + C β .
By summing up the previous inequality for i = 1, ..., N , we obtain:

1 ≤ C(m, R M ) ρ d 2 k0 e C α k 0 log 1 µ k 0 -1 N i=1 w H 1 (Bρ 0 (ζi,0)) + π L 2 (Bρ 0 (ζi,0)) + C(m, R M , N ) ρ d 2 k0 α k 0 w H 3 ( Ω) + π H 2 ( Ω) + C(N ) β ,
where Ω =

N i=1

Ωi . Moreover, µ k0-1 (d 0 + ρ 0 ) ≥ and ρ k0 = µ k0 ρ 0 imply that ρ k0 ≥ µ ρ 0 d 0 + ρ 0 .

It follows:

1 ≤ C(m, R M ) d 2 e C α k 0 log 1 µ k 0 -1 N i=1 w H 1 (Bρ 0 (ζi,0)) + π L 2 (Bρ 0 (ζi,0)) + C(m, R M , N ) w H 3 ( Ω) + π H 2 ( Ω) + C(N ) α k 0 d 2 + β . ( 31
)
Let us denote by

E = w H 3 ( Ω) + π H 2 ( Ω) + C(N ) . (32)
We simplify the last term in the right hand-side in (31) by choosing > 0 such that

α k 0 d 2 = β .
Since α k0 < 1, we obtain, for all small enough

1 ≤ C(m, R M , N ) e C α k 0 log 1 µ k 0 -1 (β+d)/α k 0 N i=1 w H 1 (Bρ 0 (ζi,0)) + π L 2 (Bρ 0 (ζi,0)) + C(m, R M , N )E β . (33)
Since 1 α k0 = e k0 log(1/α) , we obtain, using (29):

1 α k0 < e log(1/α) log((d 0 +ρ 0 )/ ) log(1/µ) +1 = 1 α d 0 + ρ 0 γ0
, where γ 0 = log(1/α) log(1/µ) . Furthermore, since 1 µ k0-1 ≤ d 0 + ρ 0 by definition of k 0 , we have

log 1 µ k0-1 < log d 0 + ρ 0 .
Then, (35) Note that the previous inequality remains true for 0 < s < s 0 , eventually by increasing the constant C(N ) involved in (32).

e c α k 0 log 1 µ k 0 -1 (β+d)/α k 0 = e 1 α k 0 c log 1 µ k 0 -1 +(β+d) log
We now look for the lower bound of this inequality with respect to s. We denote by Remembering the definition (36) of Λ, the previous inequality is equivalent to

Λ = w L 2 (Γ) + π L 2 (Γ) + ∂w ∂n L 2 (Γ) + ∂π ∂n L 2 (Γ) δ . (36
w L 2 (Γ) + π L 2 (Γ) + ∂w ∂n L 2 (Γ) + ∂π ∂n L 2 (Γ) ≥ C(M 1 , m, R M , N ).
By replacing (w, π) by ( 21), we obtain:

N j=1 |q 1 j -q 2 j | ≤ C(m, R M , M 1 , N ) u 1 -u 2 L 2 (Γ) + p 1 -p 2 L 2 (Γ) + ∂u 1 ∂n - ∂u 2 ∂n L 2 (Γ) + ∂p 1 ∂n - ∂p 2 ∂n L 2 (Γ)
which concludes the proof of Theorem 4.1, since ∂u 1 ∂n -∂u 2 ∂n = (p 1 -p 2 )n on Γ.

Figure 1 :

 1 Figure 1: Figure illustrating how informations spread in the proof of Theorem 4.1.

Figure 2 :

 2 Figure 2: Figure illustrating the construction of the sequence of balls with decreasing radius and whose center is approaching x i .

( 1 ) ≤ e 1 α

 1 ( d 0 +ρ 0 ) γ 0 (c log( d 0 +ρ 0 )+(β+d) log( 1 )) ≤ e C γ 0 log( 1 ) . (34)To summarize, for γ > γ 0 and for 0 < < 0 , we obtain from (33) and (34):1 ≤ C(m, R M , N ) e C/ γ N i=1 w H 1 (Bρ 0 (ζi,0)) + π L 2 (Bρ 0 (ζi,0)) + E β .By denoting 1 s = γ , it can be rewritten for all s > s 0 = 1 0 γ as:1 ≤ C(m, R M , N ) e Cs N i=1 w H 1 (Bρ 0 (ζi,0)) + π L 2 (Bρ 0 (ζi,0)) + E 1 s β/γ.InriaBy applying Corollary 3.5, there exist C(M 1 , m, R M , N ) > 0 and 0 < δ < 1 such that for all s > s 0 :1 ≤ C(M 1 , m, R M , N )   e Cs w L 2 (Γ) + π L 2 (Γ) + ∂w ∂n L 2 (Γ)

.

  ) and f (s) = e Cs Λ + dE 1 sβ/γ, with d ≥ 1. Let us study the function f in R * + . We have:lim s→0 f (s) = +∞, lim s→∞ f (s) = +∞.So since f is continuous on R + * , f reaches its minimum at a point s 0 > 0. At this point,f (s 0 ) = 0 ⇔ Λ = d large enough.Taking into account (37), it leads to:1 ≤ E dC(M 1 , m, R M , N ) ln E dβ CγΛ λThanks to local regularity stated in Subsection 2.2, we know that there exits a constant C(M 1 , R M ) > 0 such that:w H 3 ( Ω) + π H 2 ( Ω) ≤ C(M 1 , R M ),which, remembering the definition (32) of E, leads to E ≤ C(M 1 , R M , N ). Thus, by studying the variation of the function f y (x) = x (ln( x y )) λ on (y, +∞), for y = CγΛ dβ we obtain the existence of two positive constants C 1 (M 1 , m, R M , N ) and C 2 (M 1 , R M , N ) such that 1 ≤ C 1 (M 1 , m, R M , N ) ln C2(M1,R M ,N ) Λ λ .

RR n°8222

InriaRESEARCH CENTRE PARIS -ROCQUENCOURTDomaine de Voluceau, -Rocquencourt B.P. 105 -78153 Le Chesnay Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399


A Belkhabbaz

R Brenner

N Rupin

B Bacroix

J Fonseca

Prediction of the overall behavior of a 3D microstructure of austenitic steel by using FFT numerical scheme

The present work investigates the mechanical response of a polycrystalline aggregate by the combined use of 3D X-ray diffraction contrast tomography (DCT) experiment and the fast Fourier transform (FFT) numerical homogenisation scheme. The microstructure of a recrystallized sample of 304L austenitic steel has been characterized with the DCT technique at ESRF (ID11 beamline) which allows to have access to the spatial distribution and the individual shape of the constituent grains together with their crystalline orientation. To estimate the mechanical behaviour of such real polycrystalline microstructure, the most widely used approach is based on the finite element technique which implies a tedious and complex task of meshing. To circumvent this difficulty, an attractive approach is the FFT-based numerical scheme which allows to use the digital image of the microstructure as a direct input to compute the mechanical response of a heterogeneous material. The efficiency of this method combined with experimental tomographic data is emphasized and the elastic and elastic-plastic behavior of the polycrystalline sample is investigated. Comparisons are made with estimates delivered by the self-consistent approach, that rely on partial information on the microstructure, and first order bounds.

Introduction

Some years ago, the only available information to describe the microstructure of a metallic material was based on optical observations, to describe the shape of the grains at the surface of the material and X-ray diffraction, to describe the overall crystallographic texture of the material. During the last decade, important advances have been made with the development of automated EBSD technique [START_REF]Electron Backscatter Diffraction in Materials Science[END_REF] which gave rise to orientation imaging microscopy: the orientation and the shape of the grains at the surface of a polycrystal can be characterized. The mechanical modeling of polycrystals now makes use of these microstructural data in order to improve the description of their effective properties [START_REF] Adams | The mesostructure-properties linkage in polycrystals[END_REF] and microstructural evolutions [START_REF] Lebensohn | Orientation image-based micromechanical modelling of subgrain texture evolution in polycrystalline copper[END_REF]. However, the full-field modelling approaches generally require additional assumptions in order to obtain a 3D virtual microstructure from experimental 2D informations [START_REF] St-Pierre | 3d simulations of microstructure and comparison with experimental microstructure coming from o.i.m analysis[END_REF][START_REF] Lebensohn | Orientation image-based micromechanical modelling of subgrain texture evolution in polycrystalline copper[END_REF] and it has been shown that the morphology of the "substrate" strongly influences the fields which can be observed at the surface [START_REF] Zeghadi | Ensemble averaging stress-strain fields in polycrystalline aggregates with a constrained surface -part 2: crystal plasticity[END_REF]. A full experimental description of the 3D microstructure is thus necessary but remains difficult to obtain with the EBSD technique since it requires long and tedious treatments of the material (etching or machining) [START_REF] Rowenhorst | 3d crystallographic and morphological analysis of coarse martensite: Combining ebsd and serial sectioning[END_REF]. Even if this technique has been recently improved with the erosion process by focused ion beam (FIB) [START_REF] Matteson | Electron backscattering diffraction investigation of focused ion beam surfaces[END_REF], an inherent limitation of this method is its destructive nature. An alternative technique is the diffraction contrast tomography (DCT) which makes use of X-ray synchrotron radiation to characterize the crystalline orientation and 3D grain shape within the bulk of the polycrystal [START_REF] Johnson | X-ray diffraction contrast tomography: A novel technique for three-dimensional grain mapping of polycrystals[END_REF]. These newly experimental informations can serve to model the microstructure either with a statistical decription to enrich homogenization models or with a pointwise description to perform full field simulations. The use of these data for finite element (FE) computations, although possible [START_REF] Ludwig | New opportunities for 3d materials science of polycrystalline materials at the micrometre lengthscale by combined use of x-ray diffraction and x-ray imaging[END_REF], requires a preliminary complex task of meshing and leads to computations with a huge number of degrees of freedom. In the present work, the alternative method based on fast-Fourier transform (FFT) which has been proposed to overcome these difficulties [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] will thus be preferred. This article investigates the combined use of DCT data with the FFT numerical scheme to compute efficiently the mechanical response of a real 3D aggregate of 304L stainless steel. The elastic behavior is first studied and compared to the self-consistent estimate and first-order bounds, then some preliminary results are derived for time independent elastoplasticity.

DCT technique and 3D microstructure

The DCT [START_REF] Johnson | X-ray diffraction contrast tomography: A novel technique for three-dimensional grain mapping of polycrystals[END_REF][START_REF] Ludwig | New opportunities for 3d materials science of polycrystalline materials at the micrometre lengthscale by combined use of x-ray diffraction and x-ray imaging[END_REF] allows to capture simultaneously grains shape and orientation in the bulk of a sample. The principle of the method consists in combining 3D reconstruction from projections and diffraction techniques. The diffraction is realized thanks to a parallel, monochromatic beam of synchrotron radiation (ESRF ID11 beamline). The sample is continuously rotated through 360 • perpendicularly to the incident beam direction during the illumination. Grains orientation fulfill Bragg conditions several times during the rotation of the sample and then diffracted patterns can be observed on the screen. An algebraic reconstruction technique can be used on these indexed patterns in order to obtain each grain shape, then these constituants are reassembled in the 3D polycrystalline structure. The limitations of this technique are: it can only be used on small samples (∼500µm of diameter) of single phase material with large grains (∼20µm), the spatial resolution is of few microns and the angular precision is of 0.1 • .

The polycrystal reconstructed by DCT, composed of 362 grains, is displayed on figure 1(a) where each color is linked to a specific crystallographic orientation. For this kind of material, the mechanical anisotropy can be associated to several sources. First, each single crystal presents an anisotropic behavior which can be described by a cubic tensor in the elastic domain. The anisotropy of the polycrystalline aggregate is of course connected to the orientation of each single crystal and to the shape of the constituant. Former studies suggest that the overall anisotropy is mainly connected to the local anisotropy caracterized by Zener's parameter and to the texture of the material at least in the elastic domain. Nevertheless, the details of the mechanical fields can be affected by others microstructural features. To determine the orientation distribution function (ODF) from EBSD measurements, the harmonic method, which consists in a serie expansion of the ODF on a basis of spherical harmonic functions, has been used. The figure 1(b) shows the inverse pole figure of the austenite. It can be seen that the material exhibit a strong fiber texture in the < 111 > direction and a slightly less pronounced fiber in the < 001 >. The overall behaviour is thus expected to present a symmetry close to transverse isotropy.

Numerical homogeneization using FFT

The DCT data can be used as microstructural inputs for a micromechanical modelling scheme. As noted previously, the widely used FE technique implies a complex task of meshing for a careful description of grain boundaries and, on the opposite, the use of a simple cartesian mesh is known to have a poor numerical efficiency regarding to dedicated algorithms on this topology. Consequently, we have chosen an alternative approach to treat this problem at reduced computational cost (problems with several millions of degrees of freedom can be solved in a few minutes with this technique [START_REF] Michel | Effective properties of composite materials with periodic microstructure: a computational approach[END_REF]). Also, an attractive feature is that the image of the microstructure can be used directly as an input without using dedicated meshing algorithm. The principle is based on the fact that the local mechanical response of a periodic heterogeneous medium, discretized on a N x × N y × N z unit cell, can be computed as a integral equation. This operation convolves the Green's function of a linear reference homogeneous medium with the actual heterogeneity field, called polarization field τ τ τ(x x x). This description of the problem is similar to numerous treatments encountered in homogenization theory [START_REF] Willis | Variational and related methods for the overall properties of composites[END_REF][START_REF] Hill | Continuum micro-mechanics of elastoplastic polycrystals[END_REF] where the initial heterogeneous problem is replaced by an equivalent homogeneous problem with polarization :

σ σ σ(x x x) = C C C(x x x) : ε ε ε(x x x) becomes σ σ σ(x x x) = C C C 0 : ε ε ε(x x x) + τ τ τ(x x x) with τ τ τ(x x x) = C C C(x x x) -C C C 0 : ε ε ε (1)
If the convolution problem is difficult to solve in the real space its counterpart in the Fourier space is only a simple tensorial product. This observation has motivated the use of an algorithm based on Fast Fourier Transform (FFT) to solve this unit cell problem [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF]. In the following, the integration algorithm which follows the proposal of [START_REF] Hutchinson | Elastic-plastic behaviour of polycrystalline metals and composites[END_REF] is described hereafter. Note that the plastic slip increments ∆γ k are computed by randomly chosing a set of linearly independent and potentially active slip systems. The consistency condition delivers a relation between ∆γ k and the increment of total strain ∆ε.

Results and discussions

Young's modulus of the aggregate

In order to describe the elastic anisotropy of the overall behavior of the material, the directional dependence of Young's modulus E is computed with

1 E(n n n) = n n n ⊗ n n n : S S S : n n n ⊗ n n n (2)
Algorithm 1 Rate independent elastoplastic behaviour

σ n = g(ε n ) -strain field ε n (x g ) known, ∀ x g -elastic guess : σ E = σ n-1 + C : (ε n -ε n-1 )
-resolved shear stresses :

τ k = µ k : σ E , ∀k -plastic criterion (Schmid law) if τ k < τ k cn-1 , ∀k then σ n = σ E , ∆τ k c = 0 else -consistency condition : ∆τ k = ∆τ k c ⇒ ∆γ k , k ∈ K K = set of potentially active slip systems
-plastic strain increment :

∆ε P = k µ k ∆γ k -elastoplastic guess : σ n = σ E -C : ∆ε P -critical shear stresses increment : ∆τ k c = k h kk ∆γ k end if
where n n n is the direction of uniaxial tension. The elastic moduli characterizing the full effective response have been obtained by performing six FFT computations with independent macroscopic loadings. The effective Young's modulus obtained by FFT can be compared to results derived from analytical homogenization models, namely the self-consistent (SC) scheme and Voigt and Reuss bounds [START_REF] Willis | Variational and related methods for the overall properties of composites[END_REF]. The effective tensors for these different models are given by

C C C Voigt = C C C , ( C C C SC + C C C * ) -1 = (C C C + C C C * ) -1 , S S S Reuss = S S S . (3)
where . represents spatial average and C C C * is the "constraint tensor" introduced by Hill [START_REF] Hill | Continuum micro-mechanics of elastoplastic polycrystals[END_REF]. The results obtained are presented on figure 2. The different figures present a very similar evolution with respect to the loading direction since crystallographic texture is the first order microstructural parameter, and it is verified that the FFT and SC results is bounded by the Voigt and Reuss estimates. The obtained highest stiffness along the macroscopic z axis is consistent with the strong < 111 > fiber texture and the anisotropy of the single crystal of austenite. Interestingly, the self-consistent estimate is very close to the "exact" response obtained by FFT which takes into account the shape and spatial distribution of the individual grains in the 3D austenite sample.

Elastoplastic response and spatial distribution of the local fields

The plastic deformation occurs by glide on octahedral slip systems {111} < 110 > and the critical resolved shear stresses are constant. Figure 3 shows the uniaxial tensile response of this aggregate and the distribution of deformation ε 33 in the plastic regime. The presence of shear bands at +/-45 • spreading over many grains can be noticed. These results are in agreement with experimental observations [START_REF] Doumalin | Microextensométrie locale par corrélation d'images numériques[END_REF] and modelling results by FE [START_REF] Barbe | Three-dimensional characterization of strain localization bands in high-resolution elastoplastic polycrystals[END_REF].

Concluding remarks

In this work, DCT data have been used to compute, with a FFT numerical scheme, the effective properties and local mechanical fields of a 3D microstructure of austenitic steel. The results are found to be in good agreement with the self-consistent estimate, at least in the elastic regime, and are consistent with previous simulations and observations concerning strain localization in the plastic regime. On-going works are concerned by the implementation of more advanced constitutive relations using dislocation densities as internal variables and comparisons with other modelling approaches, namely non-linear mean field homogenization scheme and FE computations.

Figure 1 :

 1 Figure 1: Experimental characterization of the material

Figure 2 :

 2 Figure 2: Directional variation of Young's modulus for different homogenization models.

Figure 3 :

 3 Figure 3: Overall stress-strain curve of the material for unixial tensile test (on the left), spatial distribution of the ε 33 field at 1% overall strain (on the right).