Christian Licht
email: clicht@univ-montp2.fr

Somsak Orankitjaroen
email: somsak.ora@mahidol.ac.th

Patcharakorn Rojchanasuwakul

Thibaut Weller
email: weller@univ-montp2.fr

Heterogeneous linearly piezoelectric patches bonded on a linearly elastic body

Keywords:

, we studied the response of a thin homogeneous piezoelectric patch perfectly bonded to an elastic body. Here we extend this study to the case of a very thin heterogeneous patch made of a periodic distribution of piezoelectric inclusions embedded in a linearly elastic matrix and perfectly bonded to an elastic body. Through a rigorous mathematical analysis, we show that various asymptotic models arise, depending on the electromechanical loading together with the relative behavior between the thickness of the patch and the size of the piezoelectric inclusions. of B h . Let O h := ∪ S ∪ B h . The body is clamped on a part 0 of ∂ \ S with a positive two-dimensional Hausdorff measure H 2 ( 0 ) and subjected to body forces and surface forces on 1 := ∂ \ (S ∪ 0 ) of densities f and F . The whole patch is free of mechanical loading and no electric charges lay in the inclusions. We define S δ := S + δe 3 , ∀δ ∈ R and, if u s , e(u s ) and σ s denote the fields of displacement, strain and stress in O h and ϕ s , D s stand for the electric potential and the electric displacement, then part of the equations describing the electromechanical equilibriums read as

Introduction

As the essential technological interest of piezoelectric patches is the monitoring of a deformable body they are bonded to, here, in the same spirit as [START_REF] Licht | Bonding a linearly piezoelectric patch on a linearly elastic body[END_REF], we intend to propose various asymptotic models for the behavior of the body through the study of the system constituted by a very thin patch made of a periodic distribution of piezoelectric inclusions embedded in a linearly elastic matrix (studied in [START_REF] Canon | Modelling of thin elastic plates with small piezoelectric inclusions and distributed electronic circuits. Models for inclusions that are small with respect to the thickness of the plate[END_REF][START_REF] Canon | Modelling of thin isotropic elastic plates with small piezoelectric inclusions and distributed electric circuits. Models for inclusions larger or comparable to the thickness of the plate[END_REF]) perfectly bonded to a linearly elastic three-dimensional body.

A reference configuration for the linearly elastic body is a domain of R 3 assimilated to the Euclidean physical space with basis { e 1 , e 2 , e 3 }. The set lies in { x 3 < 0 } and a part of its Lipschitz-continuous boundary ∂ is a nonempty domain S of { x 3 = 0 } such that S × (-L, 0) is included in for some positive real number L. The patch perfectly bonded to the body occupies B h = S × (0, h) and is made of a periodic distribution of linearly piezoelectric inclusions perfectly bonded to a linearly elastic and electrically insulated matrix. More precisely, let Y = (0, 1) 2 and Y * be a subdomain strongly included in Y and I ε = i ∈ Z 2 ; ε(i + Y ) ⊂ S then if s = (ε, h) is the pertinent couple of geometrical parameters of the patch, the piezoelectric inclusions occupy B s I = S ε I × (0, h), S ε

I = ∪ i∈I ε ε(i + Y *
), while the elastic matrix occupies the remaining part

⎧ ⎪ ⎨ ⎪ ⎩ div σ s = f in O h u s = 0 on 0 σ s n = F on 1 σ s n = 0 on S h ∪ ∂ S × (0, h) div D s = 0 in B s I D s n = 0 on ∂ S ε I × (0, h) σ s = ae(u s ) in (σ s , D s ) = 1 h M ε (e(u s ), ∇ϕ s ) in B h (1)
f is the extension of f to B h by 0, n is the unit outward normal, a denotes the elasticity tensor, which satisfies:

a ∈ L ∞ ( ; Lin(S 3 )), ∃c > 0; c|e| 2 ≤ a(x)e • e ∀e ∈ S 3 a.e. x ∈ (2)
where Lin(S N ) is the space of linear operators on the space S N of N × N symmetric matrices where the inner product and the norm are noted • and | • | as for R 3 . If K := S 3 × R 3 is equipped with an inner product and a norm also denoted as previously and X Y * and L ∞ # (Y ; Lin(K)) denoting the characteristic function of Y * and the space of Y periodic elements of L ∞ (R 2 ; Lin(K)), respectively, then M ε is defined by:

M ε (x) = M(x/ε) ∀x ∈ B h M ∈ L ∞ # (Y ; Lin(K)) s.t. M = α -X Y * β X Y * β T X Y * γ (3) with ∃κ > 0 s.t. κ|k| 2 ≤ M( y)k • k a.e. y ∈ Y , ∀k ∈ K (4)
and, in all the sequel, for any function

w of L 2 (Y * ; R N ), X Y * w is the function defined by X Y * w = 0 in Y \ Y * , X Y * w = w in Y * .
The models indexed by p = ( p,

p 3 ), p = (p 1 , p 2 ) ∈ { 1, 2 } 2 , p 3 ∈ {
D s n = q s on S s I , (5) 1
q s being a density of charges, while p 1 = 2 corresponds to a condition of a given electrical potential:

ϕ s = ϕ s 0 on S s I (5) 2
roughly speaking p 1 = 1 deals with patches used as sensors, whereas p 1 = 2 concerns actuators (see [START_REF] Weller | Asymptotic modeling of thin piezoelectric plates[END_REF][START_REF] Reddy | On laminated composite plates with integrated sensors and actuators[END_REF]). The index p 2 accounts for the states of the interface between the piezoelectric inclusions and the body, p 2 = 1 corresponds to an electrically impermeable interface, p 2 = 2 corresponds to a grounded interface:

D s n = 0 on S ε I (6) 1 ϕ s = 0 on S ε I . (6) 2
Finally, we assume that s takes a value in a countable set with a sole cluster point such that p 3 = 1, 2, 3 if and only if lim s→0 (h/ε) = 0, 1, +∞, respectively.

It will be convenient to use the following notations

k = (ê, ĝ) ê = (e αβ ) α,β ∈{ 1,2 } ĝ = (g 1 , g 2 ) ∀k = (e, g) ∈ K k(r) = k(v, ψ) := (e(v), ∇ψ) ∀r = (v, ψ) ∈ H 1 (B s I ; R 3 ) × H 1 (B s I ) (7)
An electromechanical state will be an element r = (v, ψ) of

V p := H 1 0 (O h ; R 3 ) × p , (1,1) = H 1 m s (B s I ), (1,2) = H 1 S ε I (B s I ), (2,1) = H 1 S s I (B s I ), (2,2) = H 1 S ε I ∪S s I (B s I ) ( 8
)
where for any domain G of R

) := ϕ ∈ H 1 (B s I ); ε(i+Y * ) ϕ(x) dx = 0 ∀i ∈ I ε .
One makes the following assumptions on the data:

⎧ ⎪ ⎨ ⎪ ⎩ ϕ ε 0 denotes the restriction to S of an element of H 1/2 ({ x 3 = 0 }) still denoted by ϕ ε 0 ( f , F , q ε ) ∈ L 2 ( ; R 3 ) × L 2 ( 1 ; R 3 ) × L 2 (S), ε(i+Y * ) q ε (x) dx = 0 ∀i ∈ I ε when p = (1, 1) q s (x + he 3 ) = q ε (x), ϕ s 0 (x + he 3 ) = hϕ ε 0 (x) a.e. x ∈ S ε I (H 1 )
It is well known that for all ϕ ε

0 in H 1/2 ({ x 3 = 0 }), there exists an element of H 1 (S × (-L, 0)) when p 2 = 1, of H 1 S -L (S × (-L, 0)
) when p 2 = 2, still denoted by ϕ ε 0 , whose trace on S is ϕ ε 0 . Hence the element ϕ s 0p of p defined by

ϕ s 0p (x) = hϕ ε 0 (x, (x 3 -h)L/h) satisfies ϕ s 0p = ϕ s 0 on S s I , 1 h B s I |∇ϕ s 0p | 2 dx ≤ C (9)
Thus, classically, seeking an equilibrium state s s p leads to the problem

(P s p ) : find s s p = (u s p , ϕ s p ) in V p + (p 1 -1)(0, ϕ s 0p ) such that M s p (s s p , r) = L s p (r) ∀r ∈ V p with M s p (s, r) := ae(u) • e(v) dx + 1 h B h M ε k(s) • k(r) dx ∀s = (u, ϕ), ∀r = (v, ψ) ∈ H 1 (O h ; R 3 ) × H 1 (B s I ) (10) L s p (r) := f • v dx + 1 F • v dH 2 + (2 -p 1 ) S s I q s • ψ dx ∀r = (v, ψ) ∈ V p (11)
which, by Stampacchia's theorem, has a unique solution.

The asymptotic models

We make the following additional assumption (H 2 ) and proceed in the spirit of [START_REF] Licht | Bonding a linearly piezoelectric patch on a linearly elastic body[END_REF], but with two-scale convergence techniques instead of weak-convergence ones. Let us recall (see [START_REF] Nguetseng | A general convergence result for a functional related to the theory of homogenization[END_REF][START_REF] Allaire | Homogenization and two-scale convergence[END_REF]) that a sequence of function w ε in L 2 (S) is said to two-scale converge to a limit w 0 belonging to

L 2 (S × Y ) if, for any ψ in C ∞ 0 (S; C ∞ # (Y )), we have lim ε→0 S w ε (x)ψ(x, x/ε) dx = S Y w 0 (x, y)ψ (x, y) dx d y, ( 12
)
C ∞ # (Y )
X S ε I ϕ ε 0 two-scale converges toward X Y * ϕ 0 in L 2 (S × Y ) with ϕ 0 independent of the "fast" variable y running in Y (H 2 )
Note that this can be obtained by the strong convergence in L 2 (S) of ϕ ε 0 or by a physically realistic assumption of metallization for the upper face of the piezoelectric inclusions (see [START_REF] Canon | Modelling of thin elastic plates with small piezoelectric inclusions and distributed electronic circuits. Models for inclusions that are small with respect to the thickness of the plate[END_REF][START_REF] Canon | Modelling of thin isotropic elastic plates with small piezoelectric inclusions and distributed electric circuits. Models for inclusions larger or comparable to the thickness of the plate[END_REF]).

In the following, for any w in L 2 (B s I ; R N ) we denote its extension by zero in B h \ B s I by w.

Step 1: Convergences By arguing as in [START_REF] Licht | Bonding a linearly piezoelectric patch on a linearly elastic body[END_REF], we have:

∃c > 0; |e(u s p )| 2 dx + 1 h B h (|e(u s p )| 2 + | ∇ϕ s p | 2 ) dx ≤ C (13) which implies that k(h, s s p ) := (e(h, u s p ), g(h, ϕ s p )) := 1 h h 0 (e(u s p ), ∇ϕ s p )(•, x 3 ) dx 3 (14) satisfies |k(h, s s p )| L 2 (S;K) ≤ C (15)
thus (cf. [START_REF] Nguetseng | A general convergence result for a functional related to the theory of homogenization[END_REF][START_REF] Allaire | Homogenization and two-scale convergence[END_REF]) there exists a not relabeled subsequence such that k(h, s s p ) will two-scale converge toward some kp = (ē p , ḡp )

in L 2 (S × Y ; K) with (g p ) 3 = X Y * ϕ 0 when p = (2, 2
). In the case p 3 = 3, all the components of kp may be identified. When p 3 = 1 only some of them are clearly identified, hence, in this case, we introduce the following decomposition of K:

K = 1 K p ⊕ 2 K p ⊕ 3 K p ( 16
)
• when p = (1, 1):

1 K p = { k = (e, g) ∈ K; e i3 = 0, g 3 = 0 } , 2 K p = k = (e, g) ∈ K; ê = 0, g 3 = 0 , 3 K p = { 0 } • when p = (1, 2) or (2, 1): 1 K p = { k = (e, g) ∈ K; e i3 = 0, g = 0 } , 2 K p = k = (e, g) ∈ K; ê = 0, ĝ = 0 , 3 K p = { k = (e, g) ∈ K; e = 0, g 3 = 0 } • when p = (2, 2): 1 K p = k = (e, g) ∈ K; e i3 = 0, ĝ = 0 , 2 K p = k = (e, g) ∈ K; ê = 0, g = 0 , 3 K p = { k = (e, g) ∈ K; e = 0, g 3 = 0 }
Considering these averages in the transverse direction will be enough to deal with the cases p 3 = 2. When p 3 = 2, it suffices to proceed to a change of unknown (a so-called scaling in the mathematical theory of Kirchhoff-Love plates, see [8,[START_REF] Canon | Modelling of thin isotropic elastic plates with small piezoelectric inclusions and distributed electric circuits. Models for inclusions larger or comparable to the thickness of the plate[END_REF]) for the electromechanical state in the patch. Let h be the mapping

x = (x, x 3 ) ∈ R 3 → h x = (x, hx 3 ) ∈ R 3 (17)
which maps B 1 = S × (0, 1) onto B h , then the scaling operator is defined by:

(S h r)(( h ) -1 x) = ( h v(x), ϕ(x)) ∀r ∈ H 1 (B h ; R 3 × R) ∀x ∈ B h (18) Let s sp = (u sp , ϕ sp ) := S h s s p ,
Z = Y × (0, 1), Z * = Y * × (0, 1), z = (y, x 3 ), k z (v, ψ) := (e z (v), ∇ z ψ), 2(e z ) ij (v) := ∂ v i ∂ z j + ∂ v j ∂ z i (20) V K L = v ∈ H 1 (B 1 ; R 3 ); ∃(v M , v F ) ∈ H 1 (S; R 2 ) × H 2 (S) s.t. v(x) = v M (x) -x 3 ∇ v F (x), v 3 (x) = v F (x) (21)
By arguing as in [START_REF] Canon | Modelling of thin elastic plates with small piezoelectric inclusions and distributed electronic circuits. Models for inclusions that are small with respect to the thickness of the plate[END_REF][START_REF] Canon | Modelling of thin isotropic elastic plates with small piezoelectric inclusions and distributed electric circuits. Models for inclusions larger or comparable to the thickness of the plate[END_REF], we have Proposition 2.1.

Proposition 2.1. Let V := v ∈ H 1 0 ( ; R 3 ); γ S (v) ∈ H 1 (S; R 2 ) ,
p = γ S (u ¯p); iii) when p 3 = 3, there exists (E, G) in L 2 (S; R 3 ) × L 2 (S) and (U 1 , w 1 , 1 ) in L 2 (S; H 1 # (Y ; R 2 ) × H 1 # (Y ) × H 1 (Y * )) such that kp = k 0 p + k 1 p k 0 p = (e 0 p , g 0 p ) e 0 p = e(γ S ( ūp ), (e 0 p ) i3 = E i , 1 ≤ i ≤ 3, g 0 p = 0, (g 0 p ) 3 = G, when p = (2, 2) G = X Y * ϕ 0 k 1 p = (e 1 p , g 1 p ) e 1 p = e y (U 1 ), (e 1 p ) α3 = ∂ y α w 1 , (e 1 p ) 33 = 0, g 1 p = X Y * ∇ y 1 , (g 1 p ) 3 = 0 when p 3 = 1, there exists (U 1 , 1 ) in L 2 (S; H 1 # (Y ; R 2 ) × H 1 (Y * )) with 1 = 0 when p = (1, 1) such that 1 ( kp ) = ⎧ ⎪ ⎨ ⎪ ⎩ ( e(γ S ( ūp )) + e y (U 1 ), X Y * ∇ y 1 ) when p = (1, 1)
e( γ S ( ūp )) + e y (U 1 )

when p = (1, 2) or (2, 1)
(e( γ S ( ūp )) + e y (U 1 ), X Y * ϕ 0 ) when p = (2, 2)

when p 3 = 2, there exists (u 1 p , ϕ 1 p ) in L 2 (B 1 ; H 1 # (Y ; R 3 ) ∩ L 2 (S; H 1 (Z ; R 3 )) × L 2 (S; H 1 (Z * )) such that ê ¯p = ê(u ¯p) + êz (u ¯1 p ), (e ¯p) i3 = (e z ) i3 (û ¯p) g ¯p = X Z * ∇ z ϕ 1 p ϕ 1
p vanishes on S when p 2 = 2, ϕ 1 p = ϕ 0 on S + e 3 when p 1 = 2.

Step 2: Identification of ūp and u ¯p The process of identification of ūp and u ¯p as the unique solutions to variational problems (which implies the convergence of the whole sequences involved in Proposition 2.1) depends strongly on p 3 .

Case p 3 = 3 (i.e. h/ε → +∞): going to the limit in the variational formulation of (P s p ) with the test function

r * = (v * , ψ * ) where v * (x) = v(x) a.e. x in v(x) + min { x 3 , ε, h -x 3 } (v 1 (x, x/ε) + 2(x 3 /ε) F (x), 2w 1 (x, x/ε) + (x 3 /ε)F 3 (x)) a.e. x in B h (22) with (v 1 , w 1 ) in C ∞ 0 (S; C ∞ # (Y ; R 2 × R), F in H 1 (S; R 3
) and for all v in V the same symbol v denotes the extension into

H 1 0 (O h ; R 3 ) by a Kirchhoff-Love displacement in B h ; ψ * (x) = min { x 3 , ε, h -x 3 } (ψ 1 (x, x/ε) + (x 3 /ε)H(x)) a.e. x in B s I with ψ 1 in C ∞ 0 (S; C ∞ # (Y * )), H in H 1 (S) (23) yields ae( ūp ) • e(v) dx + S×Y M( y) kp (x, y) • (k 0 p ( γ S (v), F , H) + k 1 p (v 1 , w 1 , ψ 1 )) dx d y = f • v dx + 1 F • v dH 2 ( 24
)
where

k 0 p (v, F , H) = (e 0 p , g 0 p )(v, F , H) ∈ K s.t. e 0 p = e(v), (e 0 p ) i3 = F i , 1 ≤ i ≤ 3, g 0 p = 0, (g 0 p ) 3 = X Y * H ∀(v, F , H) in H 1 (S; R 2 ) × L 2 (S; R 3 ) × L 2 (S) k 1 p (v 1 , z 1 , ψ 1 ) = (e 1 p , g 1 p )(v 1 , z 1 , ψ 1 ) ∈ K s.t. e 1 p = e y (v 1 ), (e 1 p ) α3 = ∂ y α z 1 , (e 1 p ) 33 = 0, g 1 p = X Y * ∇ y ψ 1 , (g 1 p ) 3 = 0 ∀(v 1 , z 1 , ψ 1 ) in := H 1 # (Y ; R 2 ) × H 1 # (Y ) × H 1 (Y * )
Let i k, 1 ≤ i ≤ 7 be a basis of K eff p,3 := S 3 × {0} × {0} × R, and i k * = i k for 1 ≤ i ≤ 6 and 7 k * = X Y * 7 k, each cell problem

( i P) : find i σ in such that Y M( y)(k 1 p ( i σ ) + i k * ) • k 1 p (τ ) d y = 0 ∀τ ∈
has a unique solution (up to a constant field) and let M eff p in Lin(K eff p ) be defined by

(M eff p ) ij = Y M( y)k 1 p ( i σ + i k * ) • (k 1 p ( j σ ) + j k * ) d y 1 ≤ i, j ≤ 7 (25)
then we have:

ae( ūp ) • e(v) dx + S M eff p k 0 p ( γ S ( ūp ), E, G) • k 0 p ( v, F , H) dx = f • v dx + 1 F • v dH 2 ∀(v, F , H) ∈ H 1 (S; R 2 ) × L 2 (S; R 3 ) × L 2 (S) (26)
The additional global variables E, G may be easily eliminated through the decomposition

K eff p,3 = 1 K eff p,3 ⊕ 2 K eff p,3 (27)
with

if p = (2, 2), 1 K eff p,3 = k = (e, g) ∈ K eff p,3 ; e i3 = 0, g 3 = 0 , 2 K eff p,3 = k = (e, g) ∈ K eff p,3 ; ê = 0 (28) if p = (2, 2), 1 K eff p,3 = k = (e, g) ∈ K eff p,3 ; e i3 = 0 , 2 K eff p,3 = k = (e, g) ∈ K eff p,3 ; ê = 0, g 3 = 0 (29)
which induces a decomposition of M eff p on linear operators ij M eff p mapping i K eff p into j K eff p , and if

M eff p := 11 M eff p -12 M eff p ( 22 M eff p ) -1 21 M eff p ( 30
)
then ūp is the unique solution to

( P p,3 ) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ find ūp in V such that for all v in V ae( ūp ) • e(v) dx + S M eff p (e( γ S ( ūp ))) • e( γ S (v)) dx = f • v dx + 1 F • v dH 2 , if p = (2, 2) ae( ūp ) • e(v) dx + S M eff p (e( γ S ( ūp )), ϕ 0 e 3 ) • (e( γ S (v)), 0) dx = f • v dx + 1 F • v dH 2 , if p = (2, 2)
Case p 3 = 1 (i.e. h/ε → 0): By using successively test functions I r = ( I v, ψ), 1 ≤ I ≤ 3 and r * = (v * , ψ * ) defined by:

( I v(x)) i = x 3 δ Ii θ 1 (x)θ 2 (x/ε) ψ(x) = η(x 3 )θ 1 (x)θ 2 (x/ε) η(x 3 ) = ⎧ ⎪ ⎨ ⎪ ⎩ x 3 if p = (1, 1) or (1, 2) x 3 -h if p = (2, 1) 0 i fp = (2, 2) (θ 1 , θ 2 ) ∈ C ∞ 0 (S) × C ∞ # (Y ) δ Ii being the Kronecker symbol and v * = v in v * (x) = v(x) + εv 1 (x, x/ε) a.e. in B h where v denotes the extension into H 1 0 (O h ; R 3 ) of any element v of V by a Kirchhoff-Love displacement in B h ψ * = εψ 1 (x, x/ε) (v 1 , ψ 1 ) ∈ C ∞ 0 (S; H 1 # (Y ; R 2 )) × C ∞ 0 (S; H 1 (Y * )) we obtain ae( ūp ) • e(v) dx + S×Y M p 1 ( kp ) • 1 k(v * , ψ * ) dx d y = f • v dx + 1 F • v dH 2 (31)
where 1 k(v * , ψ * ) = ( 1 e, 1 g) is the element of 1 K p such that:

1 e = e( γ S (v)) + e y (v 1 )

1 g = ∇ y ψ 1 if p = (1, 1) ( 1 g) 3 = 0 i fp = (2, 2)
and

M p := 11 M p -12 M p ( 22 M p ) -1 21 M p ( 32
)
I J M p being the elements of the decomposition of M in linear mappings from

I K p into J K p , 1 ≤ I, J ≤ 3. Let k i , 1 ≤ i ≤ 3
be a basis of S 2 when p = (2, 2), and k i , 1 ≤ i ≤ 4 be a basis of S 2 × {0} × {0} × R when p = (2, 2), each cell problem

(P i ) : ⎧ ⎪ ⎨ ⎪ ⎩ find (w i , ψ i ) in H 1 # (Y ; R 2 ) × H 1 (Y * ) such that, for all (w, ψ) in H 1 # (Y ; R 2 ) × H 1 (Y * ) Y M p ( 1 k p (w i , ψ i ) + k i ) • 1 k p (w, ψ) d y = 0
has a unique solution (up to a constant field) and let M eff p in Lin(S 2 ) be defined by

( M eff p ) ij := Y M p ( 1 k p (w i , ψ i ) + k i ) • ( 1 k p (w j , ψ j ) + k j ) d y (33)
so that ūp is the unique solution to

( P p,1 ) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ find ūp in V such that, for all v in V ae( ūp ) • e(v) dx + S M eff p e( γ s ( ūp )) • e( γ s (v)) dx = f • v dx + 1 F • v dH 2 , if p = (2, 2) ae( ūp ) • e(v) dx + S M eff p (e( γ s ( ūp )), ϕ 0 e 3 ) • (e( γ s (v)), 0) dx = f • v dx + 1 F • v dH 2 , if p = (2, 2)
Case p 3 = 2 (i.e. h/ε → 1): as only local unknowns u 1 p , ϕ 1 p enter k ¯p , it suffices to eliminate the space variable by introducing suitable oscillating test functions in the variational formulation of (P s p ) and to go to the limit as in [START_REF] Canon | Modelling of thin isotropic elastic plates with small piezoelectric inclusions and distributed electric circuits. Models for inclusions larger or comparable to the thickness of the plate[END_REF]. Let I = {1, 2} 3 ∪ { * } with { * } = ∅ if p = (2, 2) and k I ∈ S 2 × R 3 , I ∈ I, and W be defined by

k I := ((-x 3 ) I 3 -1 e I 1 ⊗ S e I 2 , 0) if I ∈ {1, 2} 3 , k I := (0, e 3 ) if I ∈ { * } (34) W := H 1 (Z ; R 3 ) ∩ L 2 ((-1, 1); H 1 # (Y ; R 3 )) × 1 (35) 1 := ψ ∈ H 1 (Z * ); ψ = 0 on S if p 2 = 2, ψ = 0 on S + e 3 if p 1 = 2 (36)
each cell problem

(P I ) : find σ I in W such that Z M p (y)(k I + k z (σ I )) • k z (ρ) dz = 0, ∀ρ ∈ W
= f • v dx + 1 F • v dH 2 , if p = (2, 2)
Hence, we have the following convergence results

Theorem 2.1. When s goes to zero, u s p weakly converges in H 1 0 ( ; R 3 ) toward ūp , while, when p 3 = 2, u sp weakly converges in H 1 (B 1 ; R 3 ) toward u ¯p , the unique solution to (P p ).

Concluding remarks

Hence our limit models described by P p correspond to the equilibrium of the genuine body subjected to loading ( f , F ) and reinforced along S, this reinforcement being nonlocal when p 3 = 2. The characteristics of this elastic reinforcement may depend on the dielectric or piezoelectric coefficients inasmuch as these terms appear in the expression of M eff p or M eff p (see [START_REF] Weller | Asymptotic modeling of thin piezoelectric plates[END_REF], where the systematic influence of crystal symmetries has been carried out). When p = (2, 2), the applied electric potential creates a source surface term and the patch acts as an actuator. The piezoelectric inclusions being disconnected, the patch cannot act as a sensor.

  has a unique solution (up to constant fields) and let M eff p be defined by(M eff p ) I J := Z M p (y)(k I + k z (σ I )) • (k J + k z (σ J )) dz, 1 ≤ I, J ≤ (I)(37)then ( ūp , u ¯p) is the unique solution to( P p,2 ) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ find ( ūp , u ¯p) in V = ( v, v ¯) in H 1 0 ( ; R 3 ) × V K L (B 1 ) s.t γ S ( v) = (v

  3 , H 1 (G; R 3 ) and H 1 (G) respectively denote the subspaces of H 1 (G; R 3 ) and H 1 (G) of all

	elements with vanishing traces on a part	of ∂G, while H 1 m s (B s I

  ) are the subsets of H 1 (Y ) and H 1 (Y ; R 2 ) made of all elements whose traces on the opposite sides of Y are equal. So, (H 2 )

	reads as:

being the set of infinitely differentiable functions in R 2 that are periodic of period Y , while

H 1 # (Y ) and H 1 # (Y ; R 2

  we define k(h, s sp ) = (e(h, s sp ), g(h, s sp )), e(h, s sp ) = I h e(s sp )I h , I h = e 1 ⊗ e 1 + e 2 ⊗ e 2 + ⊗ e 3 , g(h, s sp ) = I h ∇s sp

	1 h	e 3 (19)
	in L 2 (B 1 ×	¯p, g ¯p)

Then k(h, s sp ) is bounded in L 2 (B 1 ; K) and, up to a subsequence, two-scale converges in L 2 (B 1 × S) to some k ¯p = (e Z ; K). In that case, one will use the additional notations:


Christian Licht

Thibaut Weller

Keywords: Asymptotic analysis, plates and rods models, piezoelectricity, patches and junctions

We aim to present mathematical models of smart devices and smart structures. Smart devices are made of materials which present significant multiphysical couplings. They are integrated in smart structures which take technological advantages of some multiphysical effects. We first propose simplified but accurate models of thin plates or slender rods made of piezoelectric or electromagneto-elastic materials in both static and dynamic cases. Then we focus on smart structures such as piezoelectric patches bonded on a linearly elastic body and piezoelectric junctions between two linearly piezoelectric or elastic bodies.

Introduction

On November 2000, Franco Maceri came to the LMGC to present the next Colloquium Lagrangianum in Taormina and a recent study [START_REF] Bisegna | A consistent theory of thin piezoelectric plates[END_REF] on piezoelectric plates. Here we present all the studies about mathematical modeling in piezoelectricity (a topic totally new for us) we did after this stimulating talk.

In the first part we intend to propose simplified but accurate models of devices made of piezoelectric or elecromagneto-elastic materials, these devices (thin plates, slender rods) presenting one or two small dimensions. We also studied the not so well-known case of piezoelectricity with electric field gradient. The models are obtained by a rigorous study of the asymptotic behavior of a three dimensional body when some of its dimensions, considered as parameters, tend to zero. We used various tools of variational and functional analysis, the point being to consider boundary value problems depending on small parameters. This study has been carried out in the steady-state and transient cases. We outline that different kinds of models appear at the limit, depending on the electrical loading. These models correspond to the physical situation when the device behaves as a sensor or as an actuator. Moreover, we are able to show that depending on the crystalline symmetry class of the material, a striking structural switch-off may appear at the limit, the device being no more piezoelectric.

The essential technological interest of piezoelectric devices being the monitoring of a deformable body they are bonded to or integrated in, the second part is devoted to smart structures. The obtained results are discussed in detail in Section 3.

Of course, this field of research has led in the past twenty years to a considerable amount of literature. In this paper, we limit ourselves to our own work. The reader will find in the references of our studies a great number of articles published in this area.

Mathematical modelings of smart devices

As usual we make no difference between the physical space and R 3 whose basis is denoted (e 1 , e 2 , e 3 ). For all ξ = (ξ 1 , ξ 2 , ξ 3 ) in R 3 , we denote (ξ 1 , ξ 2 ) by ξ. Greek indices for coordinates take their values in {1, 2} whereas latin indices run from 1 to 3.

Let H = S 3 × R 3 , where S 3 denotes the set of all 3 × 3 real and symmetric matrices. The set of all linear mappings from a space V into a space W is denoted L(V, W ) and by L(V ) if V = W .

In the sequel, for every domain G in R N , the subspace of the Sobolev space H 1 (G) whose elements vanish on Γ, included in the boundary ∂G of G, will be denoted by H 1 Γ (G).

Piezoelectric thin plates

Finding the equilibrium of a thin linearly piezoelectric plate can be formulated as follows.

The reference configuration of a linearly piezoelectric thin plate is the closure in R 3 of the set Ω ε := ω × (-ε, ε), where ω is a bounded domain of R 2 with Lipschitz boundary ∂ω and ε a small positive number. Let Γ ε lat := ∂ω × (-ε, ε), Γ ε ± := ω × {±ε} and two suitable partitions of ∂Ω ε : (Γ ε mD , Γ ε mN ) and (Γ ε eD , Γ ε eN ) with Γ ε mD and Γ ε eD of strictly positive surface measures. The plate is clamped along Γ ε mD and at an electrical potential ϕ ε 0 on Γ ε eD . It is subjected to body forces f ε in Ω ε and to surface forces g ε in Γ ε mN . Furthermore, we will consider an electrical loading d ε on Γ ε eN . We note n ε the outward unit normal to ∂Ω ε and assume that Γ ε mD = γ 0 × (-ε, ε), with γ 0 ⊂ ∂ω. The equations determining the piezoelectric state s ε := (u ε , ϕ ε ) at equilibrium are:

P(Ω ε )      div σ ε + f ε = 0 in Ω ε , σ ε n ε = g ε on Γ ε mN , u ε = 0 on Γ ε mD , div D ε = 0 in Ω ε , D ε • n ε = d ε on Γ ε eN , ϕ ε = ϕ ε 0 on Γ ε eD , (σ ε , D ε ) = M ε (x)(e(u ε ), ∇ϕ ε ) in Ω ε ,
where u ε , ϕ ε , σ ε , e(u ε ) and D ε respectively stand for the displacement, the electrical potential, the stress tensor, the tensor of small strains (i.e. the symmetrized gradient) and the electric induction. The operator M ε is an element of L(H) such that:

σ ε = a ε e(u ε ) -b ε ∇ϕ ε D ε = b ε T e(u ε ) + c ε ∇ϕ ε (1)
with b ε T the transpose of the piezoelectric tensor b ε , the elastic tensor a ε and the dielectric one c ε being symmetric and positive. Note that because of the piezoelectric coupling, M ε is not symmetric. It is easy to give a weak (or variational) formulation of the previous linear boundary problem and to conclude to the existence and the uniqueness of a solution in suitable Sobolev spaces through the Stampacchia theorem.

Nevertheless, due to the very low thickness of the plate, this classical model may be difficult to tackle numerically. The essence of our proposal of simplified but accurate modeling is to consider ε as a small parameter and to study the asymptotic behavior of s ε when ε goes to 0. In fact, two different limit behaviors indexed by p ∈ {1, 2} will occur, according to the type of boundary condition in P(Ω ε ).

From the mathematical point of view it is convenient to proceed to a change of coordinates Π ε and of unknows s p (ε) = S p (ε)s ε in order to consider functional spaces defined on a fixed domain Ω = ω × (-1, 1).

x = (x 1 , x 2 , x 3 ) ∈ Ω → Π ε x = (x 1 , x 2 , εx 3 ) ∈ Ω ε s p (ε) := (u(ε)(x), ϕ p (ε)(x)) = ((ε -1 ûε (Π ε x), u ε 3 (Π ε x)), ε -p ϕ ε (Π ε x)). (2)
The formulae defining S p (ε) stem from the assumptions on the magnitude of the electromechanical loading and are justified by the convergence results they lead to. If we consider forces and displacements, these hypotheses are the ones of [3] and supply a mathematical justification of the Kirchhoff-Love theory of thin linearly elastic plates. In addition, we assume that ϕ ε 0 has an extension into Ω ε still denoted by ϕ ε 0 and that

ϕ 0 ∈ H 1 (Ω) is such that ϕ ε 0 (Π ε x) = ε p ϕ 0 (x) with:      if p = 1 : ϕ 0 does not depend on x 3 .
if p = 2 : the closure of the projection of Γ ε eD on ω coincides with ω, moreover, either

d ε = 0 on Γ ε eN ∩ Γ ε lat or Γ ε eN ∩ Γ ε lat = ∅. (3)
Thus s(ε) is the solution of the variational problem:

Find s p (ε) ∈ (0, ϕ 0 ) + V = {r = (v, ψ) ∈ H 1 Γ mD (Ω) 3 × H 1 Γ eD (Ω)} such that Ω M (x) k p (ε, s) • k p (ε, r) dx = L(r), ∀r ∈ V
where the linear form L does not depend on ε and

     k p (ε, r) = k p (ε, (v, ψ)) = (e(ε, v), ∇ p (ε, ψ)), e(ε, v) αβ = e(v) αβ , e(ε, v) α3 = ε -1 e(v) α3 , e(ε, v) 33 = ε -2 e(v) 33 , ∇ p (ε, ψ) = ε p-1 ∇ψ, ∇ p (ε, ψ) 3 = ε p-2 ∂ 3 ψ. ( 4
)
The signs of the various powers of ε in the components of k p (ε, r) induce an orthogonal decomposition of H in subspaces H p , with ∈ {-, 0, +}, which is crucial to fully describe plates models in all admissible crystal classes. We denote by h p the projection on H p of any element h of H so that M can then be decomposed in nine elements M p ∈ L(H p , H p ), with , ∈ {-, 0, +}. Because M 00 p and M -- p are positive operators on H 0 p and H - p , the Schur complement

M p := M 00 p -M 0- p (M -- p ) -1 M -0 p (5)
is an element of L(H 0 p ). The key point of the asymptotic study is to show that if k p is the limit (in a suitable topology) of k p (ε, s p (ε)), then (M k p ) - p = (k p ) + p = 0. This will enable us to exhibit M p as the operator governing the limit constitutive equations due to the fundamental relation:

(M h) - p = h + p = 0 ⇒ M p h 0 p = (M h) 0 p and M p h 0 p • h 0 p = M h • h. (6)
The limit space of displacements will be the space of Kirchhoff-Love displacements defined by V KL := {v ∈ H 1 Γ mD (Ω) 3 ; e i3 (v) = 0} while the limit electrical spaces will be Φ e,1 := {ψ ∈ H 1 Γ eD (Ω); ∂ 3 ψ = 0} and Φ e,2 := {ψ ∈ H 1 ∂3 (Ω); ψ |Γ eD ∩Γ ± = 0}, where

H 1 ∂3 (Ω) := {ψ ∈ L 2 (Ω); ∂ 3 ψ ∈ L 2 (Ω)
} and where Γ eD stands for the image of Γ ε eD by Π ε -1 . Finally, we have the following convergence result:

Let K 1 := H 1 (Ω) and K 2 := H 1 ∂3 (Ω).
When ε → 0, the familly (s p (ε)) ε>0 of the unique solutions of P(ε, Ω) p strongly converges in X p := H 1 Γ mD (Ω) 3 × K p to the unique solution s p of

P(Ω) p

Find s ∈ (0, ϕ 0 ) + S p such that

Ω M p k(s) 0 p • k(r) 0 p dx = L(r), ∀ r ∈ S p := V KL × Φ e,p .
To get physically meaningful results, we define an electromechanical state s ε p over the real plate Ω ε by the descaling s ε p = S p (ε) -1 s p : it is the unique solution of a problem P(Ω ε ) p posed over Ω ε which is the transportation by Π ε of the (limit scaled) problem P(Ω) p . This transported problem is our proposal to model the thin linearly piezoelectric plate of thickness 2ε. Our model in fact involves two dimensional problems set on ω, which is very attractive and favourable from the numerical point of view. It is also accurate in the sense that the convergence result on the scaled states implies that s ε is asymptotically equivalent to s ε p . The first model (p = 1) with ϕ 0 = 0 deals with the physical situation when the plate is used as a sensor, the second model corresponds to an actuator. The model involves "reduced" state variables, the sole component k 0 p of the couple strain/gradient of the electrical potential, and the constitutive equation are supplied by the Schur complement (or the "condensation" of the initial operator M ε ) with respect to the maintained components. This identification is the keypoint for obtaining some decoupling and symmetry properties very important in practice (see [START_REF] Weller | Étude des symétries et modèles de plaques en piézoélectrictité linéarisée[END_REF], [START_REF] Weller | Analyse asymptotique de plaques minces linéairement piézoélectriques[END_REF] and [START_REF] Geymonat | Plates made of piezoelectric materials: When are they really piezoelectric?[END_REF]) by due account of the influence of the crystalline symmetries on the coefficients of M p . More precisely, it is possible to list some properties of the operator M p (p = 1, 2), which supplies the constitutive equations of the piezoelectric plate.

The fundamental coupling property of M remains true for M p :

M pme = -( M pem ) T , (7)
where m and e respectively denote the mechanical and electrical components of the generalized kinematics and stresses.

Considering the influence of crystalline symmetries on the three-dimensional constitutive law (see [START_REF] Royer | Ondes élastiques dans les solides[END_REF] for example), we can deduce, in the case of a polarization normal to the plate, that:

-M 2mm involves mechanical terms only, -M 1mm = M 2mm for the crystalline classes m, 32, 422, 6, 622 and 6m2, -M 1mm involves electrical terms except for these previous classes, -when p = 1, there is an electromechanical decoupling ( M pme = 0) for the classes 2, 222, 2mm, 4, 4, 422, 4mm, 42m, 6, 622, 6mm and 23, when p = 2, this decoupling occurs with the classes m, 32, 422, 6, 622 and 6m2, nevertheless the operators M pmm and M pee involve a mixture of elastic, piezoelectric and dielectric coefficients. In these cases, the plate can be considered as no more piezoelectric. We are then in a situation of a structural switch-off of the piezoelectric effect. Let us consider for example a thin piezoelectric plate constituted by a material whose crystalline symmetry class is 222. Then (1) takes the following form:

              σ 11 σ 22 σ 33 √ 2 σ 23 √ 2 σ 31 √ 2 σ 12 D 1 D 2 D 3               =               a 11 a 12 a 13 0 0 0 0 0 0 a 12 a 22 a 23 0 0 0 0 0 0 a 13 a 23 a 33 0 0 0 0 0 0 0 0 0 a 44 0 0 -b 41 0 0 0 0 0 0 a 55 0 0 -b 52 0 0 0 0 0 0 a 66 0 0 -b 63 0 0 0 b 41 0 0 c 11 0 0 0 0 0 0 b 52 0 0 c 22 0 0 0 0 0 0 b 63 0 0 c 33               •               e 11 (u) e 22 (u) e 33 (u) √ 2 e 23 (u) √ 2 e 31 (u) √ 2 e 12 (u) ϕ ,1 ϕ ,2 ϕ ,3               . ( 8
) Therefore, ( 5) leads to

M 1 =          a11a33-a 2 13 a33 a12a33-a13a23 a33 0 0 0 a12a33-a13a23 a33 a22a33-a 2 23 a33 0 0 0 0 0 a66c33+b 2 63 c33 0 0 0 0 0 c11a44+b 2 41 a44 0 0 0 0 0 c22a55+b 2 52 a5          (9)
in the sensor case and to

M 2 =      a11a33-a 2 13 a33 a12a33-a13a23 a33 0 0 a12a33-a13a23 a33 a22a33-a 2 23 a33 0 0 0 0 a 66 -b 63 0 0 b 63 c 33      (10)
in the actuator case. As outlined previsouly, the relation [START_REF] Mead | Electron transport mechanism in thin insulating films[END_REF] shows that σ and D respectively depend only on e(u) and ∇ϕ when the plate acts as a sensor, so that it can be considered as no more piezoelectric. But, when the same plate acts as an actuator, the piezoelectric coupling does not vanishes as it can be seen in [START_REF] Mindlin | Polarization gradient in elastic dielectrics[END_REF].

Electromagneto-elastic thin plates

Besides the piezoelectric coupling, some materials are sensitive to magnetic effects, thus in [START_REF] Weller | Modeling of linearly electromagneto-elastic thin plates[END_REF] we extended the previous modeling to linearly electromagneto-elastic thin plates. Now the state is described by s ε = (u ε , ϕ ε , φ ε ) where the additional variable φ ε denotes the magnetic potential and the constitutive equations read as:

σ ε = a ε e(u ε ) -b ε ∇ϕ ε -d ε ∇φ ε , D ε = b ε T e(u ε ) + c ε ∇ϕ ε + e ε ∇φ ε , B ε = d ε T e(u ε ) + e ε T ∇ϕ ε + f ε ∇φ ε . (11)
In these constitutive equations, d ε , e ε and f ε respectively stand for the piezomagnetic, electromagnetic coupling and magnetic permeability tensors, while B ε denotes the magnetic induction.

A similar mathematical analysis of the asymptotic behavior of s ε can be done to derive a simplified but accurate model of thin electromagneto-elastic plate. It involves reduced state variable and constitutive equations supplied by the condensation M ε of M ε with respect to the maintained components of (e(u ε ), ∇ϕ ε , ∇φ ε ).

But the novelty here is that four limit behaviors may appear according to the type of boundary conditions and the magnitude of the data on the electric and magnetic fields. These cases can be described as previously but by a couple of indices (p, q) ∈ {1, 2} 2 in place of the sole indice p. The physical situation when the thin plate is used as an electrical (resp. magnetic) sensor corresponds to p = 1 (resp. q = 1) while the actuator case corresponds to p, q = 2. It therefore appears two original mixed behaviors when p = q. In these situations, the plate is at the same time a sensor and an actuator excepted for the classes for which the plate is no more electromagneto-elastic (i.e. the electromechanical and magnetomechanical coefficients in M ε vanishes). The two cases p = q allow the modeling of electrically commanded magnetic devices and of magnetically commanded electric ones, which is of considerable interest in the development of non-volatile magnetic random access memories. We emphasize on the point that this behavior is here fully described for any admissible crystal class.

Piezoelectric plates with electric field gradient

In the 1960's the study of unexplained aspects of piezoelectricity led Mindlin [START_REF] Mindlin | Polarization gradient in elastic dielectrics[END_REF] to extend the classical Voigt theory [START_REF] Voigt | Lehrbuch der Kristallphysik[END_REF] in Toupin's formulation [START_REF] Toupin | The elastic dielectric[END_REF] by assuming that the stored energy function not only depends on the strain tensor and polarization vector but also on the polarization gradient tensor. What motivated Mindlin to study the effects of the polarization gradient was the capacitance of a very thin dielectric film. Experiments showed that the capacitance of a very thin film is systematically smaller than the classical prediction. Moreover, performing experimental tests, Mead [START_REF] Mead | Electron transport mechanism in thin insulating films[END_REF] showed that piezoelectric effects can also appear in centrosymmetric crystals, which is in contradiction with classical Voigt theory. And, indeed, the Mindlin's theory of elastic dielectrics with polarization gradient accomodates the observed and experimentally measured phenomena, such as electromechanical interactions in centrosymmetric materials, capacitance of thin dielectric films, surface energy of polarization, deformation and optical activity in quartz.

As in the classical piezoelectric case the physical state is described by s ε = (u ε , ϕ ε ). However, the constitutive equations read as:

σ ε = a ε e(u ε ) -b ε ∇ϕ ε -α ε ∇ 2 φ ε , D 1,ε = b ε T e(u ε ) + c ε ∇ϕ ε + β ε ∇ 2 φ ε , D 2,ε = α ε T e(u ε ) + β ε T ∇ϕ ε + γ ε ∇ 2 φ ε . (12)
Depending on the type of electric loading, three different models indexed by p appear at the limit. This result extends our previous study in [START_REF] Weller | Asymptotic modeling of thin piezoelectric plates[END_REF] and shows that gradient theory broadens the understanding of sensors and actuators. When p = 2 and p = 3 we are able to express the constitutive laws as a Schur complement of the second order piezoelectric tensor in a framework valid for any symmetry class, which means that we do not make any simplifying assumptions dealing with the crystal symmetry of the material constituting the plate. When p = 1, we are not able to explicitly derive the constitutive law of the limit model. Therefore, as in the case of first order piezoelectric rods treated treated in the next section, it seems very likely to us that non-local terms appear in this delicate situation. The study of the influence of the crystal symmetries on our models for p = 2, 3 shows that even for second order piezoelectricity, an electromechanical switch-off may appear in the structure if the plate is designed with specific materials.

Piezoelectric slender rods

From a technological point of view, piezoelectric materials can also be used in wires or slender rods. Now, the reference configuration of the piezoelectric structure is Ω ε = ε × (0, L) with L a fixed positive real number. The equations describing the equilibrium of the structure are the same as in the section 2.1 but of course the geometry of the various boundaries is different: we assume that Γ ε mD = εω × {0, L}. To get our simplified models, we proceed as in the case of plates. Due to classical assumptions on the mechanical loading (which permits the justification of Bernoulli-Navier theory of elastic slender rods (see [START_REF] Murat | Comportement asymptotique des solutions du système de l'élasticité linéarisée anisotrope hétérogène dans des cylindres minces[END_REF] and [START_REF] Trabucho | Mathematical modelling of rods[END_REF]) and on electrical loading:

if p = 1 the extension of ϕ ε 0 into Ω ε does not depend on x and Γ ε eD ⊂ εω × {0, L}, if p = 2 there exists γ e ⊂ ∂ω such that Γ ε eD ⊂ ε × (0, L). ( 13
)
the scaling is defined by:

     x = (x, x 3 ) ∈ Ω = ω × (0, L) → Π ε x = (εx, x 3 ) ∈ Ω ε s p (ε) = S p (ε)s ε (û(ε)(x), u 3 (ε)(x), ϕ(ε)(x)) = (û ε (Π ε x), ε -1 u ε 3 (Π ε x), ε -p ϕ ε (Π ε x)) (14)
so that s p (ε) is the unique solution of the variational problem

Find s p (ε) ∈ (0, ϕ 0 ) + V such that Ω M (x)k p (ε, s(ε)) • k p (ε, r) = L(r), ∀ r ∈ V,
with now:

k p (ε, (v, ψ)) = ((ε 2 e αβ (v), εe α3 (v), e 33 (v)), (ε p-2 ∇ψ, ε p-1 ∂ 3 ψ)).
As in the case of purely elastic slender rods (cf. [START_REF] Murat | Comportement asymptotique des solutions du système de l'élasticité linéarisée anisotrope hétérogène dans des cylindres minces[END_REF]), finding the limit is a little bit more difficult and the limit problems are as follows:

R(Ω) 1 Find (u, v, w, φ, ψ) ∈ V 1 such that Ω M (x) k 1 (u, v, w, φ, ψ) • k 1 (u , v , w , φ , ψ ) dx = L(u ), ∀ (u , v , w , φ , ψ ) ∈ V 1 , with                                  V 1 = V BN (Ω) × R b (Ω) × RD ⊥ 2 × Φ × Ψ, V BN (Ω) = {v ∈ H 1 Γ mD (Ω) 3 ; e αβ (v) = e α3 (v) = 0}, R b (Ω) = {v; ∃c ∈ H 1 0 (0, L); v(x) = c(x 3 )(-x 2 , x 1 ), v 3 ∈ L 2 (0, L; H 1 m (ω))}, H 1 m (ω) = {v ∈ H 1 (ω); ω ψ( x) d x = 0}, RD ⊥ 2 (Ω) = {ω; w ∈ L 2 (0, L; H 1 m (ω) 2 ), w 3 = 0 and ω (-x 2 w 1 ( x, x 3 ) + x 1 w 2 ( x, x 3 )) d x = 0, ae x 3 ∈ (0, L)}, Φ = {φ ∈ H 1 0 (0, L); φ(x) = φ(x 3 )}, Ψ = L 2 (0, L; H 1 m (ω)), k 1 (u, v, w, φ, ψ) = ( e(w), e α3 (v), e 33 (u), ∇ψ, dφ dx3 ) (15)
and

R(Ω) 2 Find (u, v, w, φ) ∈ V 2 such that Ω M (x) k 2 (u, v, w, φ) • k 2 (u , v , w , φ ) dx = L(u ), ∀ (u , v , w , φ ) ∈ V 2 ,
with

V 2 = V BN (Ω) × R b (Ω) × RD ⊥ 2 × L 2 (0, L; H 1 γe (ω)), k 2 (u, v, w, φ) = ( e(w), e α3 (v), e 33 (u), ∇φ). ( 16
)
The space V BN (Ω) is the Bernoulli-Navier displacements space.

Of course, our proposal of model is obtained by taking the inverse scaling, that is to say a transported problem R(Ω ε ) p posed over Ω ε . On the contrary to the case of plates, the state variables of the model do not reduce to the couple displacement/electrical potential but involve additional variables: two fields of displacements (easy to interpret mechanically) and a scalar field of electrical nature. Nevertheless,the kinematics of the state variables is simpler than the one of the genuine three-dimensional model which is very favourable from a numerical point of view. As in the purely elastic case it is worthwhile to note that for particular classes of monoclinic materials the additional variables v, w and ψ disappear [START_REF] Vrisysrisuwattana | Mathematical modelling of piezoelectric thin plates and slender beams through functional analysis[END_REF]. Anyway, in the case p = 1, the additional variables can be eliminated but it leads to non standard equations involving non local terms.

Dynamical response of piezoelectric plates

The interest of an efficient modeling of the dynamic response of piezoelectric plates lies in the fact that a major technological application of piezoelectric effects is the control of vibrations of structures through very thin plates or patches. We present two modelings depending on the various extents to which the magnetic effects are taken into account. Actually, because of the large discrepancy between the celerities of the mechanical and electromechanical waves, magnetic effects can be disregarded. That is why first we propose a modeling in the appropriate framework of the quasi-electrostatic approximation which claims that the electrical field still derives from an electrical potential.

Quasi-electrostatic case

Now a new parameter appears: the density ρ of the plate. In the framework of the realistic quasi-electrostatic approximation, the electrical equilibrium equation remains true but the mechanical equilibrium equation is replaced by

div σ ε + f ε = ρü ε in Ω ε
where the upper dot denotes the differentiation with respect to time. Under mild assumptions on the initial state and the essential assumption

+1 -1 x 3 M 1 (x 1 , x 2 , x 3 ) dx 3 , M 2 independent of x 3 (17)
it is possible to proceed to the study of the convergence of s ε p when ε goes to zero ( [START_REF] Weller | Étude des symétries et modèles de plaques en piézoélectrictité linéarisée[END_REF], [START_REF] Weller | Réponse dynamique asymptotique de plaques minces linéairement piézoélectriques dans l'approximation quasi-électrostatique[END_REF]), the result depends strongly on the relative behaviour of ε and ρ. A unified accurate and simplified modeling is then obtained by simply adding Ω ε ρ üε p dx to the left hand side of the equation defining the descaled limit problem P(Ω ε ) p . Thus the relationship between the reduced stress, electric displacement, strain and gradient of electrical potential remains the same as in the static case: M ε p really describes the constitutive equations of the plate! The displacement fields involved in our simplified modeling being of Kirchoff-Love type, clearly four cases, indexed by q, of relative behaviours of the parameters determine the essential nature of the limit response of the plate to the electromechanical loading:

q = 1 : ρ → ρ ∈ (0, +∞) , q = 2 : ρ → 0 and ρ/ε 2 → ∞ q = 3 : ρ/ε 2 → ρ ∈ (0, +∞) , q = 4 : ρ = o(ε 2 ). ( 18
)
In the cases q = 2 and q = 4, the limit response of the plate to the electromechanical loading is essentially quasi-static, while the cases q = 1 and q = 3 involve the acceleration of the displacement. Moreover, because of the assumption [START_REF] Voigt | Lehrbuch der Kristallphysik[END_REF], appears a decoupling between the membrane motion and the flexural one. If q = 1, 2, the flexion is neglectible and the membrane response is dynamic if q = 1, quasi-static if q = 2. When q = 3, 4, the membrane response is quasi-static whereas the flexural response is dynamic if q = 3 and quasi-static if q = 4. In these last two cases, the equation giving the flexion does not involve the limit electric potential if p = 1. The uncoupled elliptic and hyperbolic involved problems are two-dimensional and set on ω.

The steps of the derivation of our model are the following. Firts we proceed to the same scaling as in section 2.1 and to a decomposition s(ε) = s(ε) e + s(ε) r , where s(ε) e solves a problem like P(ε, Ω) p and consequently whose asymptotic behavior is provided by section 2.1. Hence s(ε) r = (u(ε) r , ϕ(ε) r ) satisfies an homogeneous variational evolution equation. Because the time derivatives do not act on ϕ(ε) r , it is possible to exhibit a linear evolution equation for u(ε) r governed by a maximal monotone operator in a suitable Hilbert space whose norm depends on (ε, ρ). Since the Trotter results of convergence of semi-goups of linear operators acting on variables spaces claim that the study of convergence of the transient problems reduces to the static case, the asymptotic behavior of u(ε) r , and consequently of s(ε) r is easily determined by straightforward variants of the convergence results of the section 2.1.

The fully dynamic case

In the previous case, the electrical field E ε was assumed to be curl-free and, consequently, equal to the gradient of the so-called electrical potential ϕ ε . If we want to take into account the magnetic effects, the state of the plate is now described by a triplet z ε = (u ε , E ε , H ε ) where H ε is the magnetic field and the equations of the problem read as:

         div σ ε + f ε = ρü ε in Ω ε Ḋε = c curl H ε in Ω ε µ Ḣε = -c curl E ε in Ω ε (σ ε , D ε ) = M ε (e(u ε ), E ε ) in Ω ε
with two kind of boundary conditions intimately linked to those of the previous cases (and, then, still indexed by p!):

p = 1 : H ε ∧ n ε = j ε on ∂Ω ε , p = 2 : E ε ∧ n ε = E ε 0 ∧ n ε on ∂Ω ε
Here c, µ, j ε , E ε 0 stand for the light celerity, the magnetic permeability, the surface current density and the exterior electrical field respectively. We will assume that there exist sufficiently smooth fields E 0 , j such that:

E ε 0 (Π ε ) = ε 2 E 0 (x), ∀x ∈ ∂Ω, j ε (Π ε ) = ε 2 j(x), ∀x ∈ Γ ± j ε α (Π ε x) = εj α (x), j ε 3 (Π ε ) = ε 2 j 3 (x), ∀x ∈ Γ lat . (19)
Let

               E ε 1 = {E ∈ L 2 (Ω ε ) 3 ; E 3 = 0, ∂ 3 E α = 0}, E ε 2 = {E; E 3 ∈ L 2 (Ω ε ); ∂ 3 (∂ α E 3 -∂ 3 E α ) = 0, E α = 0 on Γ ε ± }, H ε 1 = {H ∈ L 2 (Ω ε ) 3 ; H α = 0, ∂ 3 H 3 = 0}, H ε 2 = {H ∈ L 2 (Ω ε ) 3 ; H 3 = 0, ∂ 3 H α = 0}, Z ε p = V ε × E ε p × H ε p , k 1 (v, E) = (e αβ (v), E α ), k 2 (v, E) = (e αβ (v), E 3 ), 1 α, β 3. ( 20
) Under ( 17), [START_REF] Weller | Analyse asymptotique de plaques minces linéairement piézoélectriques[END_REF] and mild assumptions on the smoothness of the initial state, it can be shown that the state z ε p is asymptotically equivalent to

z ε p = (u ε p , E ε p , H ε p ) which satisfies:                                  Ω ε ρ üε p • v dx + Ω ε M ε p k p (u ε p , E ε p ) • k p (v, 0) dx = L ε p (v, 0), ∀v ∈ V KL (Ω ε ) ( Ḋε 1 ) 1 ( x) = -c ∂ 2 (H ε 1 ) 3 ( x) + j ε 1 ( x, ε) + j ε 1 ( x, -ε), ∀ x ∈ ω ( Ḋε 1 ) 2 ( x) = c ∂ 1 (H ε 1 ) 3 ( x) + j ε 2 ( x, ε) + j ε 2 ( x, -ε), ∀ x ∈ ω Ḋε 2 ) 2 ( x) = c (∂ 1 (H ε 2 ) 2 -∂ 2 (H ε 2 ) 1 )( x), ∀ x ∈ ω µ( Ḣε 1 ) 3 ( x) = -c (∂ 1 (E ε 1 ) 2 -∂ 2 (E ε 1 ) 1 )( x), ∀ x ∈ ω µ( Ḣε 2 ) 1 ( x) = -c (∂ 2 (E ε 2 ) 3 -∂ 3 (E ε 2 ) 2 )( x), ∀ x ∈ ω µ( Ḣε 2 ) 2 ( x) = -c (∂ 3 (E ε 2 ) 1 -∂ 1 (E ε 2 ) 3 )( x), ∀ x ∈ ω (σ ε p , D ε p ) = M ε p k p (u ε p , E ε p ) (21)
with the boundary conditions:

H ε 1 ∧ n ε = 1 2ε +ε -ε j ε (•, x 3 ) dx 3 on Γ ε lat , E ε 2 ∧ n ε = E ε 0 ∧ n ε on Γ ε ± . ( 22
)
The structure of the equations of our model is the same that those of the genuine model, but the problems are two-dimensional and with a lesser number of degrees of freedom for the state fields! Again, the key-point is to formulate a suitable scaling of the problems in terms of an evolution equation governed by a maximal monotone operator in an Hilbert space of possible states with finite scaled energy. By using Trotter theory we only have to consider the limit behavior of a perturbation of the variational equation which defines P(ε, Ω) p . This perturbation taking into account a scaling of the curl operator, the limit behavior is obtained by using weak continuity and integration by parts in the terms involving the curl operator.

Mathematical modelings of smart structures

The essential technological interest of piezoelectric devices being the monitoring of a deformable body they are bonded to or integrated in, this section is devoted to smart structures. Here we intend to propose various asymptotic models for the behavior of the body through the study of the system constituted by a very thin linearly piezoelectric flat patch perfectly bonded to or integrated in a linearly elastic or piezoelectric three dimensional body.

Piezoelectric patches

A reference configuration for the body is an open set Ω laying in { x 3 < 0 } whose part of its Lipschitz-continuous boundary ∂Ω is a non-empty domain S in { x 3 = 0 } and such that S × (-L, 0) is included in Ω for some positive real number L, while the patch occupies B ε := S × (0, ε), ε being a small real number; let O ε := Ω ∪ S ∪ B ε . The body is clamped on a part Γ 0 of ∂Ω \ S with a positive two-dimensional measure H 2 (Γ 0 ), and subjected to body forces and surface forces on Γ 1 := ∂Ω \ (S ∪ Γ 0 ) of densities f and F . Moreover, for all δ in R, let S δ denotes S + δe 3 , { e 1 , e 2 , e 3 } being a basis of the Euclidean physical space assimilated to R 3 , surface forces of density G acts on S ε whilst the patch is free of mechanical loading and electric charges in B ε and on its lateral boundary ∂S × (0, ε). If u ε , e(u ε ), σ ε denote the fields of displacement, strain and stress in O ε and ϕ ε , D ε stand for the electric potential and the electric displacement, part of the equations describing the electromechanical equilibrium read as:

         div σ ε = f in O ε , u ε = 0 on Γ 0 , σ ε n = F on Γ 1 , σ ε n = G ε on S ε , σ ε n = 0 on ∂S × (0, ε), div D ε = 0 in B ε , D ε • n = 0 on ∂S × (0, ε), σ ε = ae(u ε ) in Ω, (σ ε , D ε ) = 1 ε M e(u ε ), ∇ϕ ε in B ε , (23)
f is the extension of f to B ε by 0, n is the unit outward normal and a denotes the elasticity tensor which satisfies a ∈ L ∞ Ω; L(S 3 ) , ∃c; c|e| 2 ≤ a(x)e • e, ∀e ∈ S 3 , a.e. x ∈ Ω, (24

)
while M is an element of L ∞ S; L(H) satisfying

M = α -β β T γ , ∃κ > 0 ; κ|h| 2 ≤ M h • h, ∀h ∈ H, a.e. x ∈ S. ( 25
)
The models will be distinguished according to the additional necessary boundary conditions on S ε and S, characterized by an index p in { 1, 2 } 2 . Case p 1 = 1 corresponds to a condition for the electric displacement on S ε :

D ε • n = q ε on S ε , (26) 1
q ε being a density of electrical charges, while p 1 = 2 corresponds to a condition of given electrical potential:

ϕ ε = ϕ ε 0 on S ε ; (26) 2
roughly speaking, p 1 = 1 deals with patches used as sensors whereas p 1 = 2 concerns actuators (see [START_REF] Weller | Analyse asymptotique de plaques minces linéairement piézoélectriques[END_REF], [START_REF] Weller | Étude des symétries et modèles de plaques en piézoélectrictité linéarisée[END_REF]). Index p 2 accounts for the status of the interface between the patch and the body: p 2 = 1 corresponds to an insulating interface, p 2 = 2 corresponds to a grounded interface:

D ε • n = 0 on S, (27) 1
ϕ ε = 0 on S. (27) 2
Introducing the transverse average of the strain and of the electrical field, it is easy to go to the limit as ε goes to 0 and to show that the limit model corresponds to purely mechanical reinforcement problem along S.

Moreover, when p = (1, 2) or p = (2, 1), the electric data q or ϕ 0 does not have any influence on the limit model which corresponds to a purely elastic surface reinforcement of the body. However, the characteristics of this reinforcement may depend on the dielectric or piezoelectric coefficients (see [START_REF] Licht | Bonding a linearly piezoelectric patch on a linearly elastic body[END_REF]). On the contrary, electrical data q or ϕ 0 plays a role in models (1, 1) or (2, 2). More precisely, f , F and G being fixed, there is a one-to-one mapping between the applied electrical potential and the limit displacement. It is thus "theoretically" possible to determine what could be the electrical potential to apply on S ε in order to get a desired displacement. An approximate procedure may be done easily by finite elements. Another application is that the patch may shift the spectrum of the body in an interesting way, that is why we may regard the patch as an actuator. When p = (1, 1) there is also a one-to-one mapping between the limit displacement and the electrical charges. Thus the measurement of the latter may supply the knowledge of the state of displacements: the patch acts as a sensor!

Piezoelectric junctions

In this section, we first present our results dealing with smart structures composed of materials whose coefficients are of the same order of magnitude. However, as it is often observed that the electric permeability is very small compared to other coefficient, next we carry out a general study of piezoelectric junctions whose material coefficients are of different order of magnitude.

Piezoelectric junctions with material coefficients of same magnitude

In this section, we present various asymptotic models, indexed by p = (p 1 , p 2 ) ∈ { 1, 2, 3, 4 } 2 , for a thin piezoelectric junction between two linearly piezoelectric (p 2 = 1) or elastic (p 2 > 1) bodies. Index p 1 is relative to the magnitude of the piezoelectric coefficients of the adhesive, characterized by a single parameter µ, with respect to that of the constant thickness 2ε of a layer containing the adhesive. More precisely, we assume that h := (ε, µ) takes values in a countable set with a sole cluster point h ∈ { 0 } × [0, +∞] so that:

         p 1 = 1 : μ1 := lim h→ h(εµ) ∈ (0, +∞) p 1 = 2 : μ1 := lim h→ h(εµ) = 0, μ2 := lim h→ h(µ/2ε) = +∞ p 1 = 3 : μ2 := lim h→ h(µ/2ε) ∈ (0, +∞) p 1 = 4 : μ2 := lim h→ h(µ/2ε) = 0. ( 28
)
As previously said, index p 2 characterizes the status of the adherents but also that of the interfaces between adherents and adhesive:

        
p 2 = 1 : the two interfaces are electromechanically perfectly permeable, p 2 = 2 : the two interfaces are electrically impermeable, p 2 = 3 : one interface is electrically impermeable while the other is electroded, p 2 = 4 : the two interfaces are electroded.

(29) Let Ω be a domain, with Lipschitz-continuous boundary, whose intersection S with { x 3 = 0 } is a domain of R 2 of positive two-dimensional Hausdorff measure H 2 (S). Let Ω ± := Ω ∩ { ±x 3 > 0 } and ε be a small positive number, then adhesive and adherents occupy

B ε := S × (-ε, ε), Ω ε ± := Ω ± ± εe 3 , respectively; let Ω ε = Ω ε + ∪ Ω ε -, S ε ± := S ± εe 3 , O ε := Ω ε ∪ B ε ∪ ± S ε ± . Let (Γ mD , Γ mN ), (Γ eD , Γ eN ) be two partitions of ∂Ω with H 2 (Γ mD ), H 2 (Γ eD ) > 0 and 0 < δ := dist(Γ eD , S). For all Γ in { Γ mD , Γ mN , Γ eD , Γ eN }, Γ ± , Γ ε ± , Γ ε denotes Γ ∩ { ±x 3 > 0 }, Γ ± ± εe 3 , ∪ ± Γ ε ± , respectively; if (γ D , γ N ) is a partition of γ := ∂S, we denote { γ D , γ N , γ } × (-ε, ε) by { Γ ε DI , Γ ε NI , Γ ε lat }.
The structure made of the adhesive and the two adherents, perfectly stuck together along S ε ± , is clamped on Γ ε mD , subjected to body forces of density f ε and to surface forces of density F ε on Γ ε mN and vanishing on Γ ε lat . Moreover, a given electric potential ϕ h p0 is applied on Γ ε DI and, when p 2 = 1, on Γ ε eD , while electric charges of density d ε appear on Γ ε NI and, when p 2 = 1, on Γ ε eN . If σ h p , u h p , e(u h p ), D h p , ϕ h p stand for the fields of stress, displacement, strain, electric displacement and electric potential, respectively, the constitutive equations of the structure, for all p 1 in { 1, 2, 3, 4 }, read as:

     (σ h p , D h p ) = µM I e(u h p ), ∇ϕ h p in B ε ∀ p 2 ∈ { 1, 2, 3, 4 }, (σ h p , D h p ) = M ε E e(u h p ), ∇ϕ h p in Ω ε if p 2 = 1, σ h p = a ε E e(u h p ) in Ω ε if p 2 > 1 (30)
where

(M ε E , a ε E )(x) = (M E , a E )(x ∓ εe 3 ) ∀x ∈ Ω ε ± (31)          (M I , M E ) ∈ L ∞ S × Ω ; L(H) such that M P = a P -b P b T P c P ; ∃κ > 0, κ|k| 2 ≤ M P (x)k • k, ∀k ∈ H := S 3 × R 3 , a.e. x ∈ Ω, ∀ P ∈ { I, E } . (32)
Lastly we have to add the following conditions on S ε ± :

     p 2 = 2 D h p • e 3 = 0 on S ε ± , p 2 = 3 D h p • e 3 = 0 on S ε + , ϕ h p = ϕ h p0 on S ε -, p 2 = 4 ϕ h p = ϕ h p0 on S ε ± , (33)
the electric potential ϕ h p0 being given on S ε + or S ε ± . The same averaging method through the junction easily leads to our limit models. We are then able to show that in the case of piezoelectric adhesive and adherents (p 2 = 1), our results extend those obtained in elasticity (see [START_REF] Aitmoussa | Modélisation et études des singularités de contraintes d'un joint collé très mince[END_REF][START_REF] Licht | Some new mathematical modelings of junctions[END_REF]). The asymptotic behavior of the adhesive strongly depends on the magnitude of the stiffness compared to that of the thickness. When the magnitude of the stiffness is of the order of the inverse of the thickness, the adhesive is replaced by a material piezoelectric surface perfectly bonded to the adherents. When it is lesser, the adhesive is replaced by an electromechanical constraint between the two adherents which can be perfect adhesion, electromechanical pull-back or free separation, according to the order of magnitude of the stiffness which is, respectively, larger, equal or lower than that of the thickness.

Similarly, in the case of a thin piezoelectric layer embedded between two elastic adherents, depending on the magnitude of the stiffness, the adhesive is replaced by a material elastic surface perfectly bonded to the adherents or by a mechanical constraint between the adherents. In the case of electrically impermeable interfaces, the material surface has a non local elastic behavior (see [START_REF] Licht | Thin linearly piezoelectric junctions[END_REF]), the constitutive equations being derived from the asymptotic behavior of a thin piezoelectric plate acting as a sensor (case p = 1 in [START_REF] Weller | Asymptotic modeling of thin piezoelectric plates[END_REF]).

When one interface is electrically impermeable while the other is electroded, the material surface is an elastic membrane. When the two interfaces are electroded, the material surface is an elastic membrane with residual stress. In these last two cases, the constitutive equations are derived from the asymptotic behavior of a thin piezoelectric plate acting as an actuator (case p = 2 in [START_REF] Weller | Asymptotic modeling of thin piezoelectric plates[END_REF]). The mechanical constraint is perfect adhesion, elastic pull-back or free separation according to the order of magnitude of the stiffness. In the case of electrically impermeable interfaces, the elastic pull-back is of non local nature (since the state variable of electric nature φ, additional to the relative displacement, can be eliminated). In the two other cases, the elastic pull-back is local. When the two interfaces are electroded, it is similar to the purely elastic case, while, if only one interface is electroded, piezoelectric and dielectric coefficients enter the limit constitutive equations.

Piezoelectric hybrid junctions

Due to the wide range of values taken by the elastic, piezoelectric and dielectric coefficients of various devices, it is worthwhile to extend our previous study [START_REF] Licht | Thin linearly piezoelectric junctions[END_REF] devoted to thin linearly piezoelectric junctions to the case when the elastic, piezoelectric and dielectric coefficients of the junction are not of the same order of magnitude. Our various asymptotic models for a thin piezoelectric junction between two linearly piezoelectric or elastic bodies will be indexed by p = (p 1 , p 2 , p 3 ) in { 1, 2, 3, 4 } 3 . Indices p 1 and p 2 are respectively relative to the magnitude of the elastic and dielectric coefficients of the adhesive with respect to that of the constant thickness 2ε of the layer containing the adhesive. More precisely, we assume that h := (ε, µ) = (ε, µ mm , µ ee , µ me ) takes values in a countable set with a sole cluster point h ∈ { 0 } × [0, +∞] 3 , so that

               p 1 = 1 : μ1 mm := lim h→ h( 2εµ
               p 2 = 1 : μ1 ee := lim h→ h(2εµ ee ) ∈ (0, +∞) p 2 = 2 : μ1 ee := lim h→ h(2εµ ee ) = 0, μ2 ee := lim h→ h(µ ee /2ε) = +∞ p 2 = 3 : μ2 ee := lim h→ h(µ ee /2ε) ∈ (0, +∞) p 2 = 4 : μ2 ee := lim h→ h(µ ee /2ε) = 0. (35)
The parameters µ mm , µ ee , µ me respectively characterize the order of magnitude of the elastic, dielectric and piezoelectric coefficients of the adhesive. The case p 1 = p 2 , being already treated in [START_REF] Licht | Thin linearly piezoelectric junctions[END_REF], in the following we assume p 1 = p 2 . As in [START_REF] Licht | Thin linearly piezoelectric junctions[END_REF], index p 3 characterizes the status of the adherents but also that of the interfaces between adherents and adhesive:

        
p 3 = 1 : the two interfaces are electromechanically perfectly permeable, p 3 = 2 : the two interfaces are electrically permeable, p 3 = 3 : one interface is electrically permeable while the other one bears an electrode, p 3 = 4 : the two interfaces bear an electrode.

(36)

Therefore, the constitutive equations of the structure, for all p := (p 1 , p 2 ), read as:

     (σ h p , D h p ) = M µ I e(u h p ), ∇ϕ h p in B ε ∀p 3 ∈ { 1, 2, 3, 4 }, (σ h p , D h p ) = M ε E e(u h p ), ∇ϕ h p in Ω ε if p 3 = 1, σ h p = a ε E e(u h p ) in Ω ε if p 3 > 1 (37) where (M ε E , a ε E )(x) = (M E , a E )(x ∓ εe 3 ) ∀x ∈ Ω ε ± (38)                (M I , M E ) ∈ L ∞ S × Ω; L(H) such that M µ I := µ mm a I -µ me b I µ me b T I µ ee c I , M E := a E -b E b T E c E M P := a P -b P b T P c P ; ∃κ > 0 κ|k| 2 ≤ M P (x)k • k ∀k ∈ H a.e. x ∈ Ω, ∀P ∈ { I, E } .
(39) Lastly we have to add the following conditions on S ε ± :

     p 3 = 2 D h p • e 3 = 0 on S ε ± , p 3 = 3 D h p • e 3 = 0 on S ε + , ϕ h p = ϕ h p0 on S ε -, p = 4 ϕ h p = ϕ h p0 on S ε ± , (40)
the electric potential ϕ h p0 being given on S ε + or S ε ± . Our results show that for piezoelectric adhesive and adherents, when the elastic and dielectric coefficients of the adhesive are not of the same order, the piezoelectric coupling remains in the asymptotic model only when p = (1, 3) or (3, 1). More generally, when (necessarily only) one index p 1 or p 2 is equal to 1, the status of the limit model for the adhesive is hybrid. When p 1 = 1, the adhesive is replaced by both a material surface perfectly bonded to the adherents, from the mechanical point of view, and a constraint, from the electrical point view. On the contrary when p 2 = 1, a mechanical constraint appears with an electrical material surface perfectly permeable. The mechanical material surface is an elastic membrane with a possible nonvanishing (only when p = (1, 3)) residual stress stemming from the possible discontinuity of the electrical potential induced by the limit electrical constraint which is perfect permeability, electric pull-back or impermeability, according to the magnitude of the dielectric coefficients. The electrical material surface is of linear conductor type with a possible nonvanishing (only when p = (3, 1)) residual term stemming from the possible nonvanishing relative displacement induced by the mechanical constraint which is perfect adhesion, elastic pull-back or free separation according to the magnitude of the stiffness of the adhesive. When both p 1 and p 2 are greater than 1, the adhesive is replaced by an electromechanical constraint. As the orders of magnitude of the elastic and dielectric coefficients differ, this electromechanical constraint reduces to two independent mechanical and electrical constraints of the types previously evocated according to the values of p 1 and p 2 , respectively.

For a thin piezoelectric layer embedded between two purely elastic adherents through two electrically impermeable interfaces, the piezoelectric coupling remains in the asymptotic model only when p = (1, 3) or (3,[START_REF] Aitmoussa | Modélisation et études des singularités de contraintes d'un joint collé très mince[END_REF]. When p = (1, 3) the adhesive layer is replaced by a piezoelectric material surface; when p = (3, 1), it is replaced by a material conductive surface and a mechanical constraint. This constraint is of elastic pull-back type with a residual term stemming from the electrical potential in the conductive surface. Actually, when p 1 = 1, the adhesive layer is replaced by a material elastic surface perfectly bonded to the adherents. When p 2 = 3, the material surface has a non-local elastic behavior since the electrical potential can be eliminated, in the other cases the material surface is a standard elastic membrane. When p 1 ranges from 2 to 4, the adhesive layer is replaced by a mechanical constraint which is perfect adhesion, elastic pull-back or free separation. The elastic pull-back is nonlocal when p 2 = 1. When p 2 = 2, the electric potential vanishes, in the remaining cases the limit surface is a linear elastic conductor.

The limit models for a thin piezoelectric layer embedded between two elastic adherents, through either two electroded interfaces or one electroded and the other being impermeable, only differ when p = [START_REF] Aitmoussa | Modélisation et études des singularités de contraintes d'un joint collé très mince[END_REF]3). In all cases, there is a perfect decoupling between Electricity and Mechanics. When the magnitude of the stiffness is of the order of the inverse of the thickness, the adhesive is replaced by an elastic material membrane perfectly bonded to the adherents; when it is lesser, the adhesive is replaced by a mechanical constraint which is perfect adhesion, elastic pull-back, free separation according to the magnitude of the stiffness. The limit surface is at a given applied potential when p ∈ { 3, 4 } × { 1 }, at a vanishing one in the other cases. Actually when p = (1, 3, 3), the memory of electricity remains because piezoelectric and dielectric coefficients enter in the constitutive equations of the elastic membrane the adhesive layer reduces to.

Eventually the previous method may work when the elastic and dielectric coefficients of the junction are of the same order of magnitude with piezoelectric coefficients of lesser order. Obviously the conclusions of [START_REF] Licht | Thin linearly piezoelectric junctions[END_REF] remain but with b I replaced by 0, so that the piezoelectric coupling disappears in the asymptotic models.

  mm ) ∈ (0, +∞) p 1 = 2 : μ1 mm := lim h→ h(2εµ mm ) = 0, μ2 mm := lim h→ h(µ mm /2ε) = +∞ p 1 = 3 : μ2 mm := lim h→ h(µ mm /2ε) ∈ (0, +∞) p 1 = 4 : μ2 mm := lim h→ h(µ mm /2ε) = 0 (34)


Jérémy E Cohen

Nicolas Gillis
email: nicolas.gillis@umons.ac.be

A New Approach to Dictionary-Based Nonnegative Matrix Factorization

Keywords: dictionary, nonnegative matrix factorization, hyperspectral imaging

In this paper, we propose a new model along with an algorithm for dictionary-based nonnegative matrix factorization. We show its effectiveness on spectral unmixing of hyperspectral images using self dictionary compared to state-of-the-art methods.

I. INTRODUCTION

Low-rank matrix approximation problems have become more and more popular and have applications in a wide rang of applications; see, e.g., [START_REF] Udell | Generalized low rank models[END_REF] and the references therein. In general, the problem can be formulated as follows: given a data matrix X ∈ R p×n where each column X(:, j) is a data point in a p-dimensional space and a factorization rank r, the goal is to compute a basis matrix U ∈ R p×r and a weight matrix V ∈ R r×n such that X ≈ U V ⇐⇒ X(:, j) ≈ k k=1 U (:, k)V (k, j) ∀j.

There exists many variants of this problem, depending on the way to measure the error and the constraints on the factors (U, V ) and/or on the approximation U V .

In this paper, we study nonnegative matrix factorization (NMF) that requires that the factors U and V to be component-wise nonnegative [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF]. Moreover, we focus on a particular variant where the columns of the basis matrix U belong to a given dictionary D ∈ R p×d , where d r is the number of atoms. Mathematically, this means that U = D(:, K) for some index set K ⊂ {1, 2, . . . , d} with |K| = r. Therefore, using the Frobenius norm to quantify the error, dictionary-based NMF can be formulated as follows min K,V ≥0 ||X -D(:, K)V || 2 F such that K ⊂ {1, 2, . . . , d} and |K| = r.

An application where the model ( 1) is particularly useful for hyperspectral unmixing (HU). A hyperspectral image is an image for which there are usually between 100 and 200 channels for each pixel, corresponding to the reflectance (fraction of light reflected by that pixel) at different wavelengths. In other words, for each pixel, there is a vector containing its reflectances at different wavelengths which is its so-called spectral signature. The linear mixing model assumes that the spectral signature of each pixel is a linear combination of the spectral signatures of the constitutive materials, called endmembers, and is a valid model for macroscopic mixture of materials in the scene.

If the resolution of a hyperspectral image is high enough so that for each endmember there exists at least one pixel containing only that endmember, then the so-called pure-pixels assumption is satisfied and HU reduces to solving (1) using self dictionary D = X. A dictionary is also available when the endmembers are contained in a hyperspectral library. For more details on HU, we refer the readers to the surveys [START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches[END_REF], [START_REF] Ma | A signal processing perspective on hyperspectral unmixing: Insights from remote sensing[END_REF] and the references therein.

A. NMF with self dictionary

In this paper, we will focus on the case with self dictionary D = X. As far as we know, there are mainly two types of approaches to tackle [START_REF] Udell | Generalized low rank models[END_REF] in that case:

• Geometric approaches that selects the atoms in the dictionary based on some geometric criteria, typically based on the volume of the convex hull of X(:, K). These approaches include for example vertex component analysis (VCA) [START_REF] Nascimento | Vertex component analysis: a fast algorithm to unmix hyperspectral data[END_REF] and the successive projection algorithm (SPA) [START_REF] Araújo | The successive projections algorithm for variable selection in spectroscopic multicomponent analysis[END_REF], [START_REF] Ren | Automatic spectral target recognition in hyperspectral imagery[END_REF], [START_REF] Chan | A simplex volume maximization framework for hyperspectral endmember extraction[END_REF], [START_REF] Gillis | Fast and robust recursive algorithmsfor separable nonnegative matrix factorization[END_REF]. They are usually fast, running in O(pnr) operations. However, they do not always select atoms leading to a small data fitting term ||X -X(:, K)V || F since, most of them do not take it into account directly, as they usually put an emphasis on some geometric properties of X(:, K) (such as having a large volume). In particular, these methods are in general sensitive to outliers.

• Sparse regression approaches that are based on the following reformulation of (1), which imposes row sparsity constraints on the scores Y :

min Y ∈R d×n ||X -XY || 2 F
such that Y has r non-zero rows.

Row-sparsity of Y can be achieved in different ways; in particular using convexifications based on the 1 norm, e.g., 1,2 [START_REF] Elhamifar | See all by looking at a few: Sparse modeling for finding representative objects[END_REF], 1,∞ [START_REF] Esser | A convex model for nonnegative matrix factorization and dimensionality reduction on physical space[END_REF], or using linear programming [START_REF] Bittorf | Factoring nonnegative matrices with linear programs[END_REF], [START_REF] Gillis | Robust near-separable nonnegative matrix factorization using linear optimization[END_REF].

These methods have the advantage to better model [START_REF] Udell | Generalized low rank models[END_REF] as they take into account the data fitting term explicitly. They usually provide good solutions but are rather costly as an optimization problem in dn variables must be solved. In particular, for D = X, we have d = n hence n 2 variables. In HU, n is the order of millions and these approaches are impractical.

Hence pixels have to be selected in a preprocessing step [START_REF]A fast gradient method for nonnegative sparse regression with self dictionary[END_REF] (e.g., using a geometric approach). Moreover, the problem solved is an approximation of the original problem, which results may not be as close as desired to the solutions of the non-convex problem.

B. Contribution and outline of the paper

In this paper, we propose a new algorithm to solve the formulation (1). It will combine the advantages of the two types of approaches described above: it is fast, running in O(pnr) operations, but taking explicitly the data fitting term ||X -X(:, K)V || 2 F into account. In Section II, we describe the new algorithm, and in Section III, we show that is competes favorably with state-of-the-art approaches for HU.

II. PROPOSED ALGORITHM FOR (1)

Solving the combinatorial model (1) directly is difficult. We propose in this paper to use an alternating strategy:

• Update of V . For fixed K, solving (1) in variable V is a convex optimization problem that can be solved efficiently; it is a nonnegative least squares (NNLS) problem.

• Update of K. The update of K is difficult, being a combinatorial problem. Moreover, for V fixed, the optimal subset K is most likely to be unique hence this approach would not be able to modify an initial solution (V, K). To circumvent these difficulties, we introduce an auxiliary variable for the hidden factor U = D(:, K), and consider the problem

min K,U ≥0,V ≥0 ||X -U V || 2 F + δ||U -D(:, K)|| 2 F such that K ⊂ {1, 2, . . . , d} and |K| = r, (2)
for some penalty parameter δ > 0. Solving for U is again an NNLS problem while for K, it is now a trivial problem to pick the atoms of the dictionary the closest to the columns of matrix U . The parameter δ is progressively increased to ensure convergence of the model (2) to the model [START_REF] Udell | Generalized low rank models[END_REF]. On top of allowing to solve each subproblem effectively, the auxiliary variable U can be used as a correction of the selected atoms since U minimizes the data fitting term while being close to D(:, K) for some K. This way the self-dictionary model allows for some flexibility on the spectra actually contained in the data, which can vary from pixel to pixel due to spectral variability [START_REF] Zare | Endmember variability in hyperspectral analysis: Addressing spectral variability during spectral unmixing[END_REF].

Note that to update U and V , we do no solve the NNLS subproblems to a high precision, which is not necessary and would be computationally rather costly, but use a few steps of block coordinate descent [START_REF] Gillis | Accelerated multiplicative updates and hierarchical ALS algorithms for nonnegative matrix factorization[END_REF] (we performed 10 iterations).

Algorithm 1 provides the pseudocode for our proposed algorithm.

Algorithm 1 Alternating Optimization for (2) Input: X ∈ R p×n , integer r, maxiter. Output: U ∈ R p×r + , V ∈ R r×n +
and an index set K such that ||X -D(:, K)V || F is small and U ≈ D(:, K). 1: Choose some initial matrices U ≥ 0, V ≥ 0 and index set K, and initial value of δ. Solve for U :

min U ≥0 ||X -U V || 2 F + δ||U -D(:, K)|| 2 F .

5:

Find K so that D(:, K) matches U the best:

6: K = ∅.

7:

for k = 1 : r do 8:

K = K ∪ argmax k D(:,k) T U (:,k) ||D(:,k)||2 . 9:
end for

10: if ||U -D(:, K)|| F > 0.01||U || F then 11:
Increase δ.

12:

end if 13:
if ||U -D(:, K)|| F < 0.05||U || F and K has not changed for 5 iterations then 14:

return; the algorithm has converged.

15:

end if 16: end for a) Initialization: Our algorithm tries to solve a highly non-linear and combinatorial problem. This explains why, as we will see, it is sensitive to initialization. In this paper we use the following initialization scheme 1) Compute K using your favorite algorithm. In the numerical experiments section, we will use random initialization and several state-of-the-art pure-pixel search algorithms.

2) Define U = D(:, K) and, for V , compute the optimal solution of the unconstrained problem and project it onto the nonnegative orthant, that is, use

V = max 0, argmin Y ||X -U Y || 2 F (in Matlab, V=max(0,U\M)).
3) Improve (U, V ) using an NMF algorithm (this assumes δ = 0). We used 10 iterations of A-HALS [START_REF] Gillis | Accelerated multiplicative updates and hierarchical ALS algorithms for nonnegative matrix factorization[END_REF]. 4) Initialize

δ = 0.01 ||X -U V || 2 F ||U -D(:, K)|| 2 F
, so that the data fitting term has initially more importance in the objective function. In fact, when δ is large, the algorithm is less likely to be able to update K since U will be very close to D(:, K). b) Other related algorithms: There are a number of arbirary choices that were made to design Algorithm 1, we here shortly discuss some variants. First, in the estimation of V , D(:, K) can be swapped with U to minimize the same cost function for the two factors U and V . We noticed however that by doing so, the performance of the algorithm decreased. In particular, when D(:, K) is used in the estimation of V , after some iterations K should not change much whereas U can still be modified, so that the number of iterations is larger in this modified version. Second, the flexible approach is not mandatory. Indeed, U can be updated using a two-step procedure involving nonnegative least squares followed by a projection on the set of atoms. An advantage of this method is that it does not require to tune or fix any parameter. However we found this simpler version of the algorithm to impose too hard constraints in the first iterations where exploring the set of unconstrained U seems important. This is the reason why δ is introduced, and starts at a relatively small value.

III. NUMERICAL EXPERIMENTS ON HYPERSPECTRAL

IMAGES

In this section, we compare Algorithm 1 with three geometric algorithms, namely VCA [START_REF] Nascimento | Vertex component analysis: a fast algorithm to unmix hyperspectral data[END_REF], SPA [START_REF] Araújo | The successive projections algorithm for variable selection in spectroscopic multicomponent analysis[END_REF] and SNPA [START_REF] Gillis | Successive nonnegative projection algorithm for robust nonnegative blind source separation[END_REF], one clustering-based algorithm, H2NMF [START_REF] Gillis | Hierarchical clustering of hyperspectral images using rank-two nonnegative matrix factorization[END_REF], and a sparse regression framework from [START_REF]A fast gradient method for nonnegative sparse regression with self dictionary[END_REF], referred to as FGNSR (where a subset of 100 and 500 columns is identified using H2NMF).

We will initialize Algorithm 1 with the atoms extracted by the above methods and refer to the corresponding algorithm as d-X, where X is the algorithm; for example, d-VCA stands for Algorithm 1 initialized with VCA. We will also use 10 random initializations (picking r pixels at random as initial endmembers) and report the worst, average and best solution (in terms of reconstruction error), denoted RAND-wo, RAND-av and RANDbe, respectively. For all these numerical experiments, we increase δ using 1.5δ. In order to keep δ bounded, we only increase it if ||U -D(:, K)|| F is large; see step 10 of Algorithm 1.

When we report the CPU time of Algorithm 1, we do not report the time for the initialization. The CPU time of FGNSR does not take into account the preprocessing time by H2NMF.

We compare the different approaches on the same data sets as [START_REF]A fast gradient method for nonnegative sparse regression with self dictionary[END_REF] with the same factorization ranks:

• Urban data set with 162 wavelength and 309 × 309 pixels, r = 6, 8. • San Diego airport with 158 wavelength and 400 × 400 pixels, r = 8, 10. • The Terrain data set with 166 wavelength and 400 × 400 pixels, r = 5, 6. Tables I, II and III report the results. For an index set K extracted by an algorithm, the relative reconstruction error (rel. err.) is given by min

V ≥0 ||X -X(:, K)V || F ||X|| F .
In these tables, the relative error is given in percent. Figure 1 shows the best solution found for the Urban data set with r = 6 (with d-H2NMF), and Figure 2 the corresponding abundance maps. We can identify constitutive materials such as trees, roof tops, dirt, grass, and roads.

We observe that

• In all cases but one (FGNSR on the Urban image with r = 8, with an increase of 0.05%), Algorithm 1 is able to improve the initial solutions provided by RAND, VCA, SPA, H2NMF and FGNSR. • In all cases, Algorithm 1 converges in less than 20 iterations. The reason is that we increased δ rather aggressively.

• Even with random initial index sets K, Algorithm 1 provides solutions with small reconstruction errors.

In fact, it is rather surprising that even the worst solution generated among 10 random initial sets K is better than the solutions generated by VCA, SPA and H2NMF. Moreover, the best solution found for the San Diego airport image with r = 8 is from random initialization. This means that although Algorithm 1 is sensitive to initialization, which was expected, it allows to identify reasonable solutions regardless of the initialization. • In all cases, Algorithm 1 leads to the best solution compared to the original algorithms. A reason why Algorithm 1 works remarkably well for these data sets is because the spectral signatures of the pixels are close to one another and form a dense cloud of data points. This allows the index set K to change progressively between neighboring pixels. Analyzing the behavior of Algorithm 1 in other settings is a direction for further research.

IV. CONCLUSION

In this paper, we proposed a new algorithm for NMF with self dictionary; see Algorithm 1. Like geometric methods, it is fast, running in O(mnr) operations, hence can be applied to large problems. Like sparse-regression methods, it takes into account the data fitting term explicitly, hence identifying good atoms in the dictionary with a small approximation error. We illustrate the effectiveness of Algorithm 1 on several hyperspectral images. In all cases, it identifies the set of endmembers providing the lowest reconstruction error. We focused in this paper on this particular variant of dictionary-based low-rank matrix approximations. However, our technique can be applied to the broader class of problems, which will be presented in an extended version of this paper.

2 : 3 :

 23 for k = 1 : maxiter do Solve for V : min V ≥0 ||X -D(:, K)V || 2 F . 4:

Fig. 1 .

 1 Fig. 1. Spectral signatures of the endmembers exctracted using d-H2NMF for the Urban data set with r = 6.

Fig. 2 .

 2 Fig. 2. Abundance maps identified using d-H2NMF for the Urban data set with r = 6.

  The lowest reconstruction error is highlighted in bold. In brackets, next to the CPU time, we indicate the number of iterations needed for Algorithm 1 to converge. All tests are preformed using Matlab on a laptop Intel CORE i5-3210M CPU @2.5GHz 6GB RAM.

		r = 6		r = 8
		Time (s.)	Rel. err.	Time (s.)	Rel. err.
	RAND-wo	0.00	7.87	0.00	11.66
	d-RAND-wo	22.46 (13)	5.09	34.87 (18)	5.35
	RAND-av	0.02	11.51	0.02	9.60
	d-RAND-av	23.91 (13)	4.65	30.77 (15)	4.65
	RAND-be	0.00	13.77	0.00	5.54
	d-RAND-be	22.01 (11)	4.36	36.18 (19)	4.16
	VCA	2.01	18.38	1.86	20.11
	d-VCA	26.89 (15)	5.83	29.06 (14)	5.05
	SPA	0.30	9.58	0.30	9.45
	d-SPA	24.37 (13)	4.67	28.61 (14)	4.62
	SNPA	24.34	9.63	36.72	5.64
	d-SNPA	23.04 (13)	4.94	27.94 (13)	3.97
	H2NMF	19.02	5.81	22.35	5.47
	d-H2NMF	26.66 (15)	4.05	28.92 (14)	4.24
	FGNSR-100	2.73	5.58	2.55	4.62
	d-FGNSR-100	26.72 (14)	4.36	20.81 (8)	4.04
	FGNSR-500	40.11	5.07	39.49	4.08
	d-FGNSR-500	25.07 (13)	4.40	26.83 (12)	4.13
			TABLE I
	NUMERICAL RESULTS FOR THE URBAN DATA SET.

ACKNOWLEDGMENTS

The authors want to thank Pierre Comon, Rodrigo Cabral-Farias and Miguel A. Veganzones for helpful discussions in the early stages of this work.

The authors acknowledge the support by the F.R.S.-FNRS (incentive grant for scientific research no F.4501.16). NG also acknowledges the support by the ERC (starting grant no 679515).