Michel Weber

A Remark on Zeros of Brownian Motion

Let {W (t), t ≥ 0} be a standard Brownian motion. If I is a bounded interval on which W has no zero, an almost sure lower bound to inf{|W (t)|, t ∈ I} can be provided, when I is taken from a given countable family of intervals covering the positive half-line.

Main Result

Let {W (t), t ≥ 0} be a standard Brownian motion. Let I be some bounded interval of R + . Suppose W (t) = 0, for all t ∈ I. What can be said about the size of inf{|W (t)|, t ∈ I}? This one only depends on the location of I and of the size of I. The object of this note is to prove more precisely the following result. The proof relies upon several intermediate results on infima of |W |, which are also of independent interest.

Local infima of Brownian motion

In this section, we collect some properties of the infimum of W over bounded intervals. Precise estimates of the probability

P inf a≤t≤b |W (t) -M | ≥ c .
will be necessary. Notice preliminary, since -W and W have same law that

P inf a≤t≤b |W (t) + M | ≥ c = P inf a≤t≤b | -W (t) + M | ≥ c = P inf a≤t≤b |W (t) -M | ≥ c , so that it is enough to consider the case M ≥ 0. Put Ψ(x) = P{W (1) > x} = ∞ x e -u 2 /2 du √ 2π , x ∈ R.
The lemma below is certainly well-known, although we could not find a reference. We included a proof for the sake of completeness.

Lemma 2 Let 0 < a < b < ∞. Then for any c > 0 and any real M

P inf a≤t≤b |W (t) -M | ≥ c = |v|>c 1 -2Ψ( |v| -c √ b -a ) e -(M +v) 2 2a √ 2πa dv.
Proof. By symmetry of the law of W it suffices to consider the case M ≥ 0. By the intermediate values Theorem,

P inf a≤t≤b |W (t)-M | ≥ c = P inf a≤t≤b W (t) ≥ M +c +P sup a≤t≤b W (t) ≤ M -c . (1) Let x ≥ 0. Then P inf a≤t≤b |W (t) -M | ≥ c W (a) = M ± x = 0, if 0 ≤ x ≤ c; and if x > c, P inf a≤t≤b |W (t) -M | ≥ c W (a) = M + x = P sup a≤t≤b (W (a) -W (t)) ≤ x -c , P inf a≤t≤b |W (t) -M | ≥ c W (a) = M -x = P sup a≤t≤b (W (t) -W (a)) ≤ x -c .
As for y ≥ 0, ([1], Theorem 1.5.1)

P{ sup 0≤t≤T W (t) > y} = 2P{W (T ) > y}, we get if |x| > c, P inf a≤t≤b |W (t) -M | ≥ c W (a) = M + x = 1 -2Ψ( |x| -c √ b -a ). (2)
Therefore

P inf a≤t≤b |W (t) -M | ≥ c = R P inf a≤t≤b |W (t) -M | ≥ c W (a) = u e -u 2 2a √ 2πa du = |u-M|>c P inf a≤t≤b |W (t) -M | ≥ c W (a) = u e -u 2 2a √ 2πa du = |u-M|>c 1 -2Ψ( |u -M | -c √ b -a ) e -u 2 2a √ 2πa du = |v|>c 1 -2Ψ( |v| -c √ b -a ) e -(M +v) 2 2a √ 2πa dv, (3)
as claimed.

Remark 3 It follows from Lemma 2 that

P inf a≤t≤b |W (t) -M | > 0 = lim c↓0 P inf a≤t≤b |W (t) -M | ≥ c = R 1 -2Ψ( |v| √ b -a ) e -(M +v) 2 2a √ 2πa dv. (4) Thus P inf a≤t≤b |W (t) -M | ≥ 0 = 1 = R 1 -2Ψ( |v| √ b -a ) e -(M +v) 2 2a √ 2πa dv,
yielding a discontinuity at 0. Take for instance a = 1, b = 1 + µ 2 ; the integral above is

R 1 -2Ψ( |v| µ ) e -(M +v) 2 2 √ 2π dv → 0 µ → ∞.
And

P inf a≤t≤b |W (t) -M | = 0 = 2 R Ψ( |v| √ b -a ) e -(M +v) 2 2a √ 2πa dv.
(5)

We will also show

Lemma 4 There exists an absolute constant C, such that for every real M , and

0 < a < b < ∞, P inf a≤t≤b |W (t) -M | = 0 = 2 R Ψ( |v| √ b -a ) e -(M +v) 2 2a √ 2πa dv ≤ C min 1, b -a a e - M 2 8 max(a,b-a)
.

In particular

P inf a≤t≤b |W (t) -M | = 0 ≤ C min 1, b -a a .
Proof. By (2.6)

P inf a≤t≤b |W (t) -M | = 0 = 2 R Ψ( |v| √ b -a ) e -(M +v) 2 2a √ 2πa dv = 2 |v|≤ M 2 + |v|> M 2 Ψ( |v| √ b -a ) e -(M +v) 2 2a √ 2πa dv ≤ 2 πa e -M 2 8a R Ψ( |v| √ b -a )dv + |v|> M 2 Ψ( |v| √ b -a )dv = 2 πa e -M 2 8a √ b -a R Ψ(|w|)dw + |v|> M 2 Ψ( |v| √ b -a )dv . (6)
Recall that the Mills' ratio R

(x) = e x 2 /2 ∞ x e -t 2 /2 dt = ( √ 2π)e x 2 /2 Ψ(x) veri- fies for all x ≥ 0 R(x) ≤ π/2. Thus |v|> M 2 Ψ( |v| √ b -a )dv = √ b -a |w|> M 2 √ b-a Ψ(|w|)dw ≤ C √ b -a |w|> M 2 √ b-a e -w 2 /2 dw ≤ C √ b -ae -M 2 8(b-a) . (7)
Therefore

P inf a≤t≤b |W (t) -M | = 0 ≤ C b -a a e -M 2 8a + e -M 2 8(b-a) ≤ C b -a a e - M 2 8 max(a,b-a) . (8)
One can recover as a special case that (see [START_REF] Durrett | Stochastic calculus: A practical introduction[END_REF] p.248)

P inf a≤t≤b |W (t)| = 0 = 1 - 2 π arctan a b -a ,
or, equivalently P W (t) has no zero in (a, b) = (2/π) arcsin a/b. It is possible to also give an exact expression of the probability P inf a≤t≤b |W (t) -M | = 0 , although for M = 0 this one is relatively more complicated. This is indicated in the Lemma below.

Lemma 5 Let 0 < a ≤ b < ∞. We have

P inf a≤t≤b |W (t)-M | = 0 = - 2 π M √ a 1 √ b-a b ue -(M u) 2 2a |x|≤M 1-u 2 a e -x 2 2 dx √ 2π du + 1 - 2 π e -M 2 2a arctan a b -a .
In particular

P inf a≤t≤b |W (t)| = 0 = 1 -2 π arctan a/(b -a).
And for every positive real c

P 0 < inf a≤t≤b |W (t)| < c = 2 c/ √ a 0 1 -2Ψ( u √ a √ b -a ) e -u 2 2 √ 2π du +4 ∞ c/ √ a Ψ( u √ a -c √ b -a ) -Ψ( u √ a √ b -a ) e -u 2 2 √ 2π du. Proof. Let M 1 = M/ √ a, and put F (s) = R 2Ψ |w|s e -(M 1 +w) 2 2 dw √
2π . We have

P inf a≤t≤b |W (t) -M | = 0 = R 2Ψ( |v| √ b -a ) e -(M +v) 2 2a √ 2πa dv = R 2Ψ |w| a b -a e -(M 1 +w) 2 2 √ 2π dw = F a b -a . ( 9
)
As ∂ ∂s Ψ(|w|s) = |w|Ψ ′ (|w|s) = -|w| √ 2π e -(|w|s) 2 /2 , we have ∂ ∂s F (s) = R 2 ∂ ∂s Ψ (|w|s e -(M 1 +w) 2 2 √ 2π dw = - 2 π R |w|e -1 2 [(ws) 2 +(M1+w) 2 ] dw √ 2π . But [(ws) 2 + (M 1 + w) 2 ] = [w √ s 2 + 1 + M1 √ s 2 +1 ] 2 + M 2 1 s 2 s 2 +1 , hence ∂ ∂s F (s) = - 2 π e - M 2 1 s 2 2(s 2 +1) R |w|e -1 2 [w √ s 2 +1+ M 1 √ s 2 +1 ] 2 dw √ 2π = - 2 π e - M 2 1 s 2 2(s 2 +1) s 2 + 1 R |z|e -1 2 [z+ M 1 √ s 2 +1 ] 2 dz √ 2π = - 2 π e - M 2 1 s 2 2(s 2 +1) s 2 + 1 E g - M 1 √ s 2 + 1 , (10)
where g denotes a Gaussian standard random variable. But for any real a

E |g + a| = |a| |a| -|a| e -x 2 /2 dx √ 2π + 2 π e -a 2 /2 . ( 11
)
Hence

∂ ∂s F (s) = - 2 π e - M 2 1 s 2 2(s 2 +1) s 2 + 1 M 1 √ s 2 + 1 |x|≤ M 1 √ s 2 +1 e -x 2 /2 dx √ 2π + 2 π e - M 2 1 2(s 2 +1) = - 2 π M 1 e - M 2 1 s 2 2(s 2 +1) (s 2 + 1) 3/2 |x|≤ M 1 √ s 2 +1 e -x 2 /2 dx √ 2π - 2 π e -M 2 1 /2 s 2 + 1 . (12)
If M = 0, this takes a much simplified form

∂ ∂s F (s) = - 2 π 1 s 2 + 1 . Further F (0) = R e -(M 1 +w) 2 2 dw √ 2π = 1. Therefore F (s) -1 = - 2 π M 1 s 0 e - M 2 1 t 2 2(t 2 +1) |x|≤ M 1 √ t 2 +1 e -x 2 /2 dx √ 2π dt (t 2 + 1) 3/2 - 2 π e -M 2 1 /2 s 0 dt t 2 + 1 = -e -M 2 1 /2 2 π M 1 s 0 |x|≤ M 1 √ t 2 +1 e -x 2 /2 dx √ 2π e M 2 1 2(t 2 +1) dt (t 2 + 1) 3/2 - 2 π arctan s . (13)
Consequently,

P inf a≤t≤b |W (t) -M | = 0 = F a b -a = - 2 π M √ a 1 √ b-a b ue -(M u) 2 2a |x|≤M 1-u 2 a e -x 2 2 dx √ 2π du + 1 - 2 π e -M 2 2a arctan a b -a . ( 14
)
If M = 0, this is simplified into

P inf a≤t≤b |W (t)| = 0 = 1 -2 π arctan a b-a . ( 15
)
As concerning the second formula, it suffices to write

P 0 < inf a≤t≤b |W (t)| ≤ c = P inf a≤t≤b |W (t)| ≤ c -P inf a≤t≤b |W (t)| = 0 = 1 -2 ∞ c √ a 1 -2Ψ( u √ a -c √ b -a ) e -u 2 2 √ 2π du -1 -2 ∞ 0 1 -2Ψ( u √ a √ b -a ) e -u 2 2 √ 2π du = -2 ∞ c √ a 1 -2Ψ( u √ a -c √ b -a ) e -u 2 2 √ 2π du +2 ∞ 0 1 -2Ψ( u √ a √ b -a ) e -u 2 2 √ 2π du = 2 c √ a 0 1 -2Ψ( u √ a √ b -a ) e -u 2 2 √ 2π du +4 ∞ c √ a Ψ( u √ a -c √ b -a ) -Ψ( u √ a √ b -a ) e -u 2 2 du √ 2π . Notice that π 2 = ∞ 0 dt (1+t 2 )
, and so

1 - 2 π arctan s = 2 π ∞ 0 dt (1 + t 2 ) - s 0 dt (1 + t 2 ) = 2 π ∞ s dt (1 + t 2 ) ≤ C min(1, s -1 ).
We deduce the bound obtained in Lemma 2.2

P inf a≤t≤b |W (t)| = 0 ≤ C min 1, b -a a .
In what follows we shall be interested in finding estimates of the delicate random variable

β -α [a,b] • χ {β [a,b] >0} , 0 ≤ α < 1, where we set β [a,b] := inf a≤t≤b |W (t)|.
Proposition 6 Let b > a > 0, and η > 0. Then,

P 0 < β [a,b] ≤ η ≤ 16η √ 2π min 1 √ b -a , 1 √ a .
Further, for any real α, 0 ≤ α < 1

E 1 β α [a,b] • χ{β [a,b] > 0} ≤ 47 1 -α min 1 √ b -a , 1 √ a + 1.
The result gives a control which is uniform in b-a. We have β . By using the second formula in Lemma 5

P 0 < β ≤ η = 4 η/ √ a 0 1 -2Ψ( u √ a √ b -a ) e -u 2 2 √ 2π du +4 ∞ η/ √ a Ψ( u √ a -η √ b -a ) -Ψ( u √ a √ b -a ) e -u 2 2 √ 2π du. (16) As 1 -2Ψ(x) = x -x e -u 2 /2 du √ 2π ≤ min((2/π) 1/2 x, 1), x ≥ 0 η/ √ a 0 1 -2Ψ( u √ a √ b -a ) e -u 2 2 √ 2π du ≤ η/ √ a 0 min 2 π 1/2 u √ a √ b -a , 1 e -u 2 2 √ 2π du ≤ min 2 π 1/2 √ a √ b -a η/ √ a 0 u e -u 2 2 √ 2π du , η √ 2πa ≤ min η max u≥0 ue -u 2 2 π √ b -a , η √ 2πa = η min 1 π √ b -a , 1 √ 2πa . (17) Furthermore Ψ( u √ a -η √ b -a ) -Ψ( u √ a √ b -a ) = u √ a √ b-a u √ a-η √ b-a e -u 2 2 du √ 2π ≤ η 2π(b -a) , (18)
which implies

∞ η/ √ a Ψ( u √ a -η √ b -a ) -Ψ( u √ a √ b -a ) e -u 2 2 √ 2π du ≤ η 2π(b -a) ∞ η/ √ a e -u 2 2 √ 2π du ≤ η 2π(b -a) . (19)
Besides, with the variable change u

= v (b-a) a , letting v η = η √ b-a . ∞ η/ √ a Ψ( u √ a -η √ b -a ) -Ψ( u √ a √ b -a ) e -u 2 2 √ 2π du = b -a a ∞ vη Ψ(v -v η ) -Ψ(v) e -v 2 (b-a) 2a √ 2π dv.
We have b -a a

2vη vη Ψ(v -v η ) -Ψ(v) e -v 2 (b-a) 2a √ 2π dv ≤ b -a a v η √ 2π = η √ 2πa . Now if v ≥ 2v η , then v -v η ≥ v -v/2 = v/2. Consequently b -a a ∞ 2vη Ψ(v -v η ) -Ψ(v) e -v 2 (b-a) 2a √ 2π dv = b -a a ∞ 2vη v v-vη e -x 2 2 √ 2π dx e -v 2 (b-a) 2a √ 2π dv ≤ b -a a v η ∞ 2vη e -v 2 /8 √ 2π e -v 2 (b-a) 2a √ 2π dv ≤ b -a a v η ∞ 0 e -v 2 /8 dv 2π = η √ 2πa . It follows that ∞ η/ √ a Ψ( u √ a -η √ b -a ) -Ψ( u √ a √ b -a ) e -u 2 2 √ 2π du ≤ min 2η √ 2πa , η √ 2π(b-a) ≤ 2η √ 2π min 1 √ a , 1 √ b-a . ( 20
)
By reporting

P 0 < β ≤ η ≤ 4η min 1 π √ b-a , 1 √ 2πa + 8η √ 2π min 1 √ a , 1 √ b-a ≤ 16η √ 2π min 1 √ b-a , 1 √ a . (21)
Now let X be a random variable such that X ≥ 0 a.s. As χ{]0

, 1]} = ∞ n=0 χ{] 1 2 n+1 , 1 2 n ]}, we have by the Beppo-Levi Theorem E 1 X α •χ{X > 0} = ∞ n=0 E 1 X α •χ{ 1 2 n+1 < X ≤ 1 2 n } +E 1 X α •χ{X > 1}, ( 22
)
where "=" means that both expresions are simultaneously finite or infinite. And

∞ n=0 E 1 X α • χ{ 1 2 n+1 < X ≤ 1 2 n } + E 1 X α • χ{X > 1} ≤ ∞ n=0 2 α(n+1) P 0 < X ≤ 1 2 n + 1. ( 23
)
Apply this with X = β. Then E β -α • χ{β > 0} will be finite once we prove that the series ∞ n=0 2 α(n+1) P 0 < β ≤ 1 2 n is convergent. By using (2.13) with η = 2 -n , we get

∞ n=0 2 α(n+1) P 0 < β ≤ 1 2 n ≤ 2 α+4 √ 2π min 1 √ b -a , 1 √ a ∞ n=0 2 -n(1-α) ≤ 32 2 1-α -1 min 1 √ b -a , 1 √ a ≤ 47 (1 -α) min 1 √ b -a , 1 √ a . ( 24
)
And we conclude that

E 1 β α • χ{β > 0} ≤ 47 1 -α min 1 √ b -a , 1 √ a + 1,
as claimed.

3 Proof of Theorem

Theorem 1

 1 Let ϑ k ≥ 0 be such that T N = k≤N ϑ k ↑ ∞ and denoteI N = [T N , T N +1 ]. Let η k ≥ 0 be such that )| ≥ η N or W (t) = 0 f or some t ∈ I N , N ultimately = 1.

  [a,b] → 0 as b-a → ∞, but in the same time the constraint β [a,b] > 0 becomes also stronger, making the probability of the set {β [a,b] > 0} small. When b -a → 0, β [a,b] → |W (a)|, and this is reflected by the term 1/ √ a in our estimate. Proof. Write β = β [a,b]

  Recall thatI N = [T N , T N +1 ]. It is now easy. As N ≥1 P 0 < β IN ≤ η N ≤ C Borel-Cantelli Lemma that P inf t∈IN |W (t)| ≥ η N or W (t) = 0 for some t ∈ I N , N ultimately = 1.

Acknowledgments. I thank Istvan Berkes for valuable remarks and careful reading of the paper.


Patcharakorn Viriyasrisuwattana

Christian Licht

Somsak Orankitjaroen

Thibaut Weller

Thin hybrid linearly piezoelectric junctions Les jonctions minces hybrides linéairement piézoélectriques

Keywords:

 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.

r é s u m é

Nous étendons notre étude [1] consacrée aux jonctions minces linéairement piézoélec-triques au cas où les coefficients élastiques, piézoélectriques et diélectriques de la jonction ne sont pas du même ordre de grandeur.

Introduction

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εμ 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 ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩
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

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

The physical situation is that of [START_REF] Licht | Thin linearly piezoelectric junctions[END_REF], which we recall as follows. Let be a domain, with Lipschitz-continuous boundary, of R 3 , assimilated with the physical Euclidean space with basis { e 1 , e 2 , e 3 }, 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 , eD ), ( 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 and subjected to body forces of density f ε and to surface forces of density F ε on ε mD that vanishes on ε lat . Moreover, a given electric potential ϕ h p 0 is applied on ε DI (and also on ε eD when p 3 = 1), while electric charges of density d ε appear on ε NI (and also on ε

eN when p 3 = 1). If σ h p , u h p , e(u h p ), D h p , ϕ h
p respectively stand for the fields of stress, displacement, strain, electric displacement and electric potential, 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 (3)
where

(M ε E , a ε E )(x) = (M E , a E )(x ∓ εe 3 ) ∀x ∈ ε ± (4) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ (M I , M E ) ∈ L ∞ S × ; Lin(K) 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 ∈ K := S 3 × R 3 a.e. x ∈ , ∀P ∈ { I, E } (5)
and Lin(K) is the space of linear operators on K whose inner product and norm are noted • and | • | as in R 3 (the same notations for the norm and inner product also stand for S N the space of N × N symmetric matrices).

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 p 0 on S ε - p 3 = 4 ϕ h p = ϕ h p 0 on S ε ± ( 6
)
the electric potential ϕ h p 0 being given on S ε + or S ε ± . 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 ∈ H 1 (O; R 3 × R) e(v) ∈ D (S; S 2 ); e(v) αβ = 1 2 (∂ α v β + ∂ β v α ), 1 ≤ α, β ≤ 2, ∀ v ∈ D (S; R 3 ) (7)
and the same symbol e(•) shall also stand for the symmetrized gradient in the sense of distributions of D (O; R 3 ), O ∈ O ε , , \ S, B ε , ε or D (S; R 2 ). An electromechanical state with vanishing electric potential on ε DI and on ε eD when

p 3 = 1 will belong to V ε p := H 1 ε mD (O ε ; R 3 ) × ε p 3 , with ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ε 1 = H 1 ε DI ∪ ε eD (O ε ) ε 2 = H 1 ε DI (B ε ) if H 2 ( ε DI ) > 0, H 1 m (B ε ) if H 2 ( ε DI ) = 0 ε 3 = H 1 ε DI ∪S ε - (B ε ) ε 4 = H 1 ε DI ∪ ± S ε ± (B ε ) (8)
where, for any domain

O of R N , N = 2, 3, H 1 (O; R M ) denotes the subspace of H 1 (O; R M ), M = 1
⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ Given ( f , F , d E , d I ) in L 2 ( ; R 3 ) × L 2 ( mN ; R 3 ) × L 2 ( eN ) × L 2 γ N × (-1, 1) with lat d I dH 2 = 0 when p 3 = 2 and H 2 γ D × (-1, 1) = 0 ϕ oI in H 3/2 (R) vanishing in { |x 3 | > 1 + δ/2 }, and ϕ oE in H 1 ( ) vanishing on S, then: f ε (x) = f (x ∓ εe 3 ) a.e. x ∈ ε ± , f ε (x) = 0 a.e. x ∈ B ε F ε (x) = F (x ∓ εe 3 ) a.e. x ∈ ε mN± d ε (x) = (2μ ee ) 1/2 d I (x, x 3 /ε) a.e. x ∈ ε NI d ε (x) = d E (x ∓ εe 3 ) a.e. x ∈ ε eN± if p 3 = 1 ϕ h p o (x) = ϕ oE (x ∓ εe 3 ) + ε p DI ϕ oI (x ± (1 -ε)e 3 ) a.e. x ∈ ε ± ε p DI ϕ oI (x, x 3 /ε) a.e. x ∈ B ε (9)
where p DI is such that

p DI = 0 if ∂ 3 ϕ oI = 0 in S × (-1, 1), p DI = 1 if ∂ 3 ϕ oI = 0 in S × (-1, 1). We also introduce the element ϕ o of H 1,1 ( , S) := ψ ∈ H 1 ( ) whose trace γ o (ψ) on S belongs to H 1 (S) defined by ϕ o (x) = ϕ oE (x) + (1 -p DI )ϕ oI (x ± e 3 ) a.e. x ∈ ± . We note ϕ o the trace on γ D of ϕ o and set ϕ oI = 1 2 ϕ oI (•, 1) -ϕ oI (•, -1) .
Then, if M p and L p are defined by:

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ M p (s, r) := ε M ε E k(s) • k(r) dx + B ε M μ I k(s) • k(r) dx, if p 3 = 1 ε a ε E e(u) • e(v) dx + B ε M μ I k(s) • k(r) dx, if p 3 > 1 L p (r) := f ε • v dx + ε mN F ε • v dH 2 + ε N I ∪ ε d ε ψ dH 2 ε = ε eN if p 3 = 1, ε = ∅ if p 3 > 1 ( 10
)
seeking an equilibrium state leads to the problem

(P h p ) : Find s h p in (0, ϕ h p 0 ) + V ε p such that M p (s h p , r) = L p (r), ∀ r ∈ V ε p
which, by Stampacchia's theorem, has a unique solution.

The asymptotic models

By proceeding as in [START_REF] Licht | Thin linearly piezoelectric junctions[END_REF], we will determine the asymptotic behavior of the structure when h goes to h under the following assumption on the behavior of μ me , whose rationale will clearly appear in Step 3 below.

(H p ) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ There exists μme in [0, +∞) such that μme = lim h→ h μ me , with lim h→ h μ 2 me μ ee 1 ε = 0 when p = (4, 1, 4) or p = (3, 2, p 3 ), (3, 4, p 3 ), (4, 2, p 3 ), 1 ≤ p 3 ≤ 4 lim h→ h μ 2 me μ ee ε = 0 when p = (2, 1, 4) or p = (1, 2, p 3 ), (1, 4, p 3 ), (2, 4, p 3 ), 1 ≤ p 3 ≤ 4 lim h→ h μ 2 me μ mm 1 ε = 0 when p = (2, 3, p 3 ), (2, 4, p 3 ), (4, 3, p 3 ), 1 ≤ p 3 ≤ 4 lim h→ h μ 2 me μ mm ε = 0 when p = (4, 2, 1) or p = (2, 1, p 3 ), (4, 1, p 3 ), 1 ≤ p 3 ≤ 4 lim h→ h μ me ε = 0 whenp = (4, 3, 4) lim h→ h μ me = 0 whenp = (1, 3, 4) or p = (2, 3, 4) lim h→ h μ me = μme when p = (1, 3, p 3 ), 1 ≤ p 3 ≤ 3, or p = (3, 1, p 3 ), 1 ≤ p 3 ≤ 4
In the following, C will denote various constants independent of h which may vary from line to line. It will be convenient in the cases p 3 > 1 to use the same symbol s h p for (u h p , φh p ) where φh p denotes the extension into ε of ϕ h p by 0. Without loss of generality, we suppose H 2 ( mD+ ) > 0; moreover, we assume H 2 ( mD-) > 0 when p 1 = 4, and H 2 ( eD± ) > 0 when p 2 = 4.

Step 1 (a priori estimates): By taking r = s h p -(0, ϕ h p 0 ) in the variational formulation of (P h p ), one has:

μ mm |e(u h p )| 2 L 2 (B ε ;S 3 ) + μ ee |∇ϕ h p | 2 L 2 (B ε ;R 3 ) + |k(s h p )| 2 L 2 ( ε ;K) ≤ C (11)
Step 2 (convergence of (s h p )): As in [START_REF] Licht | Thin linearly piezoelectric junctions[END_REF], the two following tools are suitable to describe the asymptotic behavior of the electromechanical state in the adherents and adhesive, respectively. First, let T ε be the mapping from

H 1 ( ε ; R 3 × R) into H 1 ( \ S; R 3 × R) defined by: (T ε r)(x) = T ε (v, ψ) (x) = (T ε 1 v, T ε 2 ψ)(x) := (v, ψ)(x ± εe 3 ) ∀x ∈ ± (12) Note that T ε s h p = (T 1 u h p , 0) if p 3 > 1! For any w in H 1 ( \ S; R N ), N ∈ { 1, 3 }, if γ ± o (w ± ) denotes the trace on S of its restriction w ± to ± , JwK stands for γ + o (w + ) -γ - o (w -). Next for all r = (v, ψ) in H 1 (B ε ; R 3 × R), we set the following element of L 2 (S; K): k p (ε, r) = e p (ε, v), g p (ε, ψ) := 1 (2ε) q 1 ε -ε e(v)(•, x 3 ) dx 3 , 1 (2ε) q 2 ε -ε ∇ψ(•, x 3 ) dx 3 q i = max(2 -p i , 0), i = 1, 2 (13)
and there holds

k p (ε, r) = (e( U p ), ∇ p ), (U h p , h p ) := 1 (2ε) q 1 ε -ε u h p (•, x 3 ) dx 3 , 1 (2ε) q 2 ε -ε ϕ h p (•, x 3 ) dx 3 (14)
So (11) and standard estimates in Sobolev spaces (see [START_REF] Licht | Thin linearly piezoelectric junctions[END_REF]) imply:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ |k(T ε s h p )| L 2 ( \S;K) ≤ C , |JT ε s h p K| 2 L 2 (S;R 3 ×R) ≤ C ε 1 + 1 μ mm + 1 μ ee |e p (ε, u h p )| 2 L 2 (S;S 3 ) ≤ C ε -2q 1 ε μ mm , |g p (ε, ϕ h p )| 2 L 2 (S;R 3 ) ≤ C ε -2q 2 ε μ ee |U h p | 2 L 2 (S;R 3 ) ≤ C ε 2(1-q 1 ) 1 + ε μ mm ε 2(q 2 -1) | h p | 2 L 2 (S) ≤ Cc * (h), c * (h) = 1 + ε μ ee if p 3 = 1, 1 εμ ee if p 3 = 2, ε 2 + ε μ ee if p 3 > 2 |e( U h p )| 2 L 2 (S;S 2 ) ≤ C 1 ε 2q 1 • ε μ mm , |∇ h p | 2 L 2 (S;R 2 ) ≤ C 1 ε 2q 2 • ε μ ee |U h p -γ ± o ((T ε 1 u h p ) ± )| 2 L 2 (S;R 3 ) ≤ C ε + ε μ mm if p 1 = 1 | h p -γ ± o ((T ε 2 ϕ h p ) ± )| 2 L 2 (S) ≤ C ε + ε μ ee if p 2 = 1 (15)
Thus, if a ⊗ S b denotes the symmetrized tensor product of a and b in R 3 , we deduce: Proposition 2.1.

1.

There exists sp = ( ūp , φp ii) when p 2 = 2, p = 0 and ḡp = 0; iii) when p 2 = 3, p and ḡp vanish only when p 3 = 2, while ( ḡp ) 3 = Jϕ p K when p 3 = 1, ( ḡp ) 3 = 0 when p 3 = 3.

) in H 1 mD ( \ S; R 3 ) × H 1 eD ( \ S) such that T ε s h p weakly converges in H 1 ( \ S; R 3 × R) toward some sp = ( ū, φp ) ; φp = 0 when p 3 > 1,
As in the next step, we will show that ( ūp , φp ) is the unique solution of a variational problem, note that the whole sequences converge. When both ēp and ḡp are defined we set kp = (ē p , ḡp ).

Step 3 (identification of (s p , ēp , ḡp )): We proceed in two different ways depending on whether kp does not exist, kp is fully or partially identified.

When p = (3, 1, 4) or 1 < p 1 , p 2 ≤ 4, p 1 = p 2 , 1 ≤ p 3 ≤ 4, it suffices to go to the limit in the variational formulation of (P h p ) by using the following test functions r ε p = (v ε p , ψ ε p ) and taking duly account of the estimates (15), Proposition 2.1, Cauchy Schwarz inequality and (H p ) which are constructed in order that lim

h→ h B ε b I ∇ϕ h p • e(v ε p ) dx = lim h→ h B ε b T I e(u h p ) • ∇ψ ε p dx = 0
The test functions r ε p reads as:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ v ε p = w min (p 1 -1,2),ε , 1 ≤ p 3 ≤ 4 ψ ε p = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ ζ min (p 2 -1,2),ε p 3 = 1 (θ 1 + x 3 θ 2 )/ε, θ 1 , θ 2 ∈ C ∞ 0 (S) p 2 = 3, 4, p 3 = 2 (1 + x 3 /ε)θ, θ ∈ C ∞ 0 (S) p 2 = 3, 4, p 3 = 3 0 i f (2 ≤ p 2 ≤ 4, p 3 = 4) or (p 2 = 2, p 3 = 2, 3)
where for all

w 1 in H 1 mD ( ; R 3 ) and all ζ 1 in H 1 eD ( ) vanishing in a neighborhood of γ D , let (w 1,ε , ζ 1,ε ) be defined by w 1,ε , ζ 1,ε (x) = (w 1 , ζ 1 )(x ∓ εe 3 ) a.e. x ∈ ε ± (w 1 , ζ 1 )(x, 0) a.e. x ∈ B ε
For all w 2 in H 1 mD ( \ S; R 3 ) and all ζ 2 in H 1 eD ( \ S) vanishing in a neighborhood of γ D , let (w 2,ε , ζ 2,ε ) be defined by

(w 2,ε , ζ 2,ε )(x) = (w 2 , ζ 2 )(x ∓ εe 3 ) a.e. x ∈ ε ± (w a , ζ a )(x, x 3 /ε) + |x 3 | ε (w s , ζ s )(x, x 3 /ε) a.e. x ∈ B ε with (w a , ζ a )(x) = 1 2 [(w 2 , ζ 2 )(x, x 3 ) -(w 2 , ζ 2 )(x, -x 3 )] (w s , ζ s )(x) = 1 2 [(w 2 , ζ 2 )(x, x 3 ) + (w 2 , ζ 2 )(x, -x 3 )] When p = (1, 3) or (3, 1), with 1 ≤ p 3 ≤ 3, let M Ip := μL 1 mm a I -μme b I μme b T I μL 2 ee c I , L i = 1 + p i /2 , i = 1, 2 (16)
with • the floor function. We first prove

(M Ip kp ) 2 p = 0 ( 17
)
where k i p denotes the projection on K i p of any element k of K with

K = K 1 p ⊕ K 2 p ⊕ K 3 p (18)
K 1 p being made of the elements k = (e, g) whose nonvanishing components e ij , g l are the nonvanishing components (ē p ) ij , ( ḡp ) l of kp which are identified by Proposition 2.1 in terms of ūp , p or ϕ o , K 2 p is made of the elements whose nonvanishing components are the components of kp that are not identified by Proposition 2.1, and K 3 p is made of the elements whose nonvanishing components are the vanishing components of kp identified by Proposition 2.1. By using suitable test functions ρ ε p , we may deduce (see [START_REF] Weller | Étude des symétries et modèles de plaques en piézoélectricité linéarisée[END_REF][START_REF] Licht | Asymptotic modeling of thin piezoelectric plates[END_REF]):

(M Ip kp ) 1 = M Ip ( kp ) 1 ; M Ip := (M Ip ) 11 -(M Ip ) 11 (M Ip ) 12 (M Ip ) 22 -1 (M Ip ) 21 (19) with (M Ip ) ij , 1 ≤ i, j ≤ 3,
ρ ε p (x) = (x 3 + ε)(I p 1 w, I p 2 p 3 ψ(x)) a.e. x ∈ B ε 2ε(I p 1 w + , I p 2 p 3 ψ + )(x -εe 3 ) in ε + , 0 in ε - ( 20
)
where, given

(w, ψ) in C ∞ o (S; R 3 × R), (w + , ψ + ) is an extension into H 1 mD+ ( + ; R 3 ) × H 1 eD+ ( + ), and I p 1 = max(0, 2 -p 1 ), I 1p 3 = 1 if p 3 ≤ 2, I 1p 3 = 0 if p 3 > 2, I 3p 3 = 0 if p 3 = 1, I 3p 3 = 1 if p 3 > 1. Second, given (v, ψ) in H 1 mD ( ; R 3 ) × H 1 eD ( ) ∩ H 2 ( ; R 3 × R), ψ vanishing in a neighborhood of γ D , we define r ε p = (v ε p , ψ ε p ) by:
when p 1 = 1:

v ε p (x) = v(x, 0) -x 3 ∇ v 3 (x, 0), v ε 3 (x) = v 3 (x, 0) a.e. x ∈ B ε v ε p (x) = v(x ∓ εe 3 ) ∓ εR ± ∇ v 3 (•, 0), 0 (x ∓ εe 3 ) a.e. x ∈ ε ± , when p 1 = 3: v ε p = w 2,ε when p 2 = 1: ⎧ ⎪ ⎨ ⎪ ⎩ ψ ε p (x) = ψ(x ∓ εe 3 ) in ε ± , ψ(x, 0) in B ε if p 3 = 1 ψ ε p (x) = ψ(x, 0) in B ε if p 3 = 2 ψ ε p (x) = 0 in B ε if p 3 ≥ 3 , when p 2 = 3 : ψ ε p = ⎧ ⎪ ⎨ ⎪ ⎩ ζ 2,ε if p 3 = 1 ψ(x, 0) if p 3 = 2 0 i fp 3 = 3 (21)
where R ± is a continuous lifting operator from

H 1/2 (S; R 3 ) into H 1 mD± ( ± ; R 3 ). Hence (H p ) and (19) yield: lim h→ h B ε M μ I k(s h p ) • k(r ε p ) dx = S M Ip ( kp ) 1 • (e , g ) dx (22) where:
when p = (1, 3) : e = e( v), g = Jζ Ke

3 if p 3 = 1, g = ∇ψ if p 3 = 2, g = 0 if p 3 = 3 when p = (3, 1) : e = Jw 2 K ⊗ S e 3 , g = ∇ψ for 1 ≤ p 3 ≤ 3 (23)
In the remaining cases, as, respectively, ḡp = 0 or does not exist, or ēp = 0 or does not exist, we proceed in the same way, but with a suitable decomposition of S 3 or R 3 , respectively, and M Ip replaced by μ1 mm ãI or μ1 ee cI , respectively, ãI and cI being defined in a similar way as M Ip .

Lastly, Jensen inequality and the previously established weak convergences achieve the proof of the following convergence result, which supports our asymptotic models in the form of variational problems P p .

Theorem 2.1.

• If p 3 = 1, when h goes to h, T ε s h p strongly converges in H 1 ( \ S; R 3 × R) toward sp the unique solution to P ( p,1) :

Find s = (u, ϕ) in (0, ϕ o ) + V p 1 × p 2 1 such that M ( p,1) (s, r) = L ( p,1) (r) ∀r = (v, ψ) ∈ V p 1 × p 21
where

M ( p,1) (s, r) = M E k(s) • k(r) dx + M I( p,1) (s, r) M I( p,1) (s, r) := ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 0 if p = (2, 4) or (4, 2) S μ1 mm ãI e( û) • e( v) dx if p = (1, 2) or (1, 4) S μ1 ee cI ∇ϕ • ∇ψ dx if p = (2, 1) or (4, 1) S μ2 ee c I JϕKe 3 • JψKe 3 dx if p = (2, 3) or (1, 3) S μ2 mm a I JuK ⊗ S e 3 • JvK ⊗ S e 3 dx if p = (3, 2) or (3, 4) S M Ip (e( û), JϕKe 3 ) • (e( v), JψKe 3 ) dx if p = (1, 3) S M Ip JuK ⊗ S e 3 , ∇γ o (ϕ) • JvK ⊗ S e 3 , ∇γ o (ψ) dx if p = (3, 1) L p (r) := ⎧ ⎨ ⎩ f • v dx + mN F • v dH 2 + eN d E ψ dH 2 + ( μ1 ee ) 1/2 γ N 1 -1 d I (•, x 3 ) dx 3 ψ dl if p 2 = 1 f • v dx + mN F • v dH 2 + eN d E ψ dH 2 if p 2 ≥ 2 V 1 := v ∈ H 1 mD ( ; R 3 ); v ∈ H 1 (S; R 2 ) , V 2 := H 1 mD ( ; R 3 ), V 3 = V 4 := H 1 mD ( \ S; R 3 ) 11 := ψ ∈ H 1 eD ( ); γ o (ψ) ∈ H 1 γ D (S
Find s = (u, φ) in 0, q 2 γ o (ϕ o ) + V p 1 × p 2 p 3 such that M p s, r = L p (r) ∀r = (v, ψ) ∈ V p 1 × p 2 p 3
where M p s, r := a e(u) • e(v) dx + M Ip (s, r)

• p 3 = 2 M Ip (s, r) := ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 0 if p = (2
• p 3 = 3 M Ip (s, r) := ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 0 if p ∈ { 2, 4 } × {
• p 3 = 4 M Ip (s, r) := ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 0 if p ∈ { 2, 4 } × { 1, 3, 4 } S M Ip (JuK ⊗ S e 3 , ∇γ o (ϕ o ), ϕ oI ) • (JvK ⊗ S e 3 , 0) dx if p = (3, 1) S μ2 mm a I JuK ⊗ S e 3 • JvK ⊗ S e 3 dx i f p = (3, 2) or (3, 4) S μ1 mm ãI e( û) • e( v) dx i f p ∈ { 1 } × { 2, 3, 4 } L p (r) := ⎧ ⎨ ⎩ f • v dx + mN F • v dH 2 + ( μ1 ee ) 1/2 γ N 1 -1 d I (•, x 3 ) dx 3 ψ dl if p 2 = 1 f • v dx + mN F • v dH 2 if p 2 ≥ 2 (p 2 ,2) := H 1 γ D (S) or H 1 m (S) according to the positivity of the length of γ D (p 2 ,3) = (p 2 ,4) := { 0 } , p 2 = 2 (2,p 3 ) := { 0 } , 2 ≤ p 3 ≤ 4

Concluding remarks

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 [START_REF] Licht | Asymptotic modeling of thin piezoelectric plates[END_REF][START_REF] Licht | Thin linearly piezoelectric junctions[END_REF]. 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 [START_REF] Licht | Asymptotic modeling of thin piezoelectric plates[END_REF][START_REF] Licht | Thin linearly piezoelectric junctions[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's 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 = (1, 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 potential ϕ o 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 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 piezoelectric coupling disappears in the asymptotic models.

  and ūp belongs to H 1 ( ; R 3 ) when p 1 ≤ 2, φp belongs to H 1 ( ) when p 2 ≤ 2 and p 3 = 1. 2. When p 1 = 4, e p (ε, u h when p 2 = 1, p is equal to γ o ( φp ) when p 3 = 1 or to γ o (ϕ o ) when p 3 ≥ 3, furthermore the trace on γ D of p is equal to φo while ḡp = ∇ p and ( ḡp ) 3 = ϕ oI when p 3 = 4;

p ) weakly converges in L 2 (S; S 3 ) toward some ēp , and there exists U p in H 1 (S; R 3 ) such that U h p weakly converges in H 1 (S; R 2 ) toward U p , (U h p ) 3 strongly converges in L 2 (S) toward (U p ) 3 , moreover i) when p 1 = 1, U p = γ 0 ( ūp ), ēp = e( ūp ); ii) when p 1 > 1, U p = 0, and ēp = Ju p K ⊗ S e 3 . 3. When p 2 = 4, g p (ε, ϕ h p ) weakly converges in L 2 (S; R 3 ) toward some ḡp and there exists p in H 1 (S) such that h p weakly converges in H 1 (S) toward p ; moreover, i)

  being the decomposition of M Ip in linear operators mapping K i p into K

	using (H p ) and ρ ε	j p . That is obtained by

p defined by:

  ( \ S; R 3 ) toward ūp while h p converges, strongly in H 1 (S) if p 2 ≤ 3 and strongly in L 2 (S) if p 3 = 4 and lim h→ h ε 3 /μ ee = 0, toward p the unique solution to P

	) ,	21 := H 1	eD ( ),	31 = 34 := H 1	eD ( \ S)
	• If p 3 > 1, when h goes to h, T ε 1 u h p strongly converges in H 1

p :

  , 4) or (4, 2) JuK ⊗ S e 3 • JvK ⊗ S e 3 dx if p = (3, 2) Ip JuK ⊗ S e 3 , ∇φ • JvK ⊗ S e 3 , ∇ψ dx if p = (3, 1)

	S mm a I S μ2 μ2 ee c I ∇φ • ∇ψ dx if p = (2, 3)

S μ1 mm ãI e( û) • e( v) dx if p = (1, 2) or (1, 4) S μ1 ee cI ∇φ • ∇ψ dx if p = (2, 1) or (4, 1)

S M Ip (e( û), ∇φ) • (e( v), ∇ψ) dx if p = (1, 3)

S M

  Ip (JuK ⊗ S e 3 , ∇γ o (ϕ o ), ϕ oI ) • (JvK ⊗ S e 3 , 0) dx if p = (3, 1)

			1, 3, 4 }
	S	μ2 mm a I JuK ⊗ S e 3 • JvK ⊗ S e 3 dx	i f p = (3, 2) or (3, 4)
	S	μ1 mm ãI e( û) • e( v) dx	i f p = (1, 2) or (1, 4)
	S M Ip e( û) • e( v) dx	i f p = (1, 3)

S M


P Taforel

M Renouf

F Dubois

J.-C Voivret

FEM/DEM coupling strategies for the modelling of ballast/soil interaction

Keywords: Discrete Element Method (DEM), Finite Element Method (FEM), hybrid FEM/DEM model, Non-Smooth Contact Dynamics (NSCD), 3D railway track modelling, Soil/ballast interaction, LMGC90

Numerous modelling strategies have been used over the last decade to improve the understanding of the physical behaviour of track system thanks to numerical simulation. In this context, Discrete Element Methods (DEM) are traditionally used to describe the granular behaviour of the ballast layer of the railway track whereas Finite Element Methods (FEM) are used to model the whole track system (soil+ballast) in which ballast is modelled as a continuous media. If both approaches offer relevant modelling frames to study railway track, this works aims to go further by coupling those two approaches in an hybrid FEM/DEM model making possible to improve our understanding of the global behaviour of the railway track.

Introduction

Looking for performance and reliability and at the same time controlling computational costs are the main issues in the railway community. Industrials, scientists or infrastructure managers have indeed to deal with many aspects surrounding track systems from which it results a technological innovation race so as to propose improved services to users while insuring their security.

In this context, improving the understanding of the behaviour of the track system has never been more important than today. Thanks to the development of high performance modelling tools, numerical simulations are commonly used in this way. As an essential component of the track system, the ballast -which plays a major role in the railway track regarding the stress distribution, the vibration damping, the lateral anchorage of sleepers, etc -is particularly studied. Much more than understanding the behaviour of this coarse granular material itself, this is its behaviour in the whole track system -in other words its interaction with the other components of the track (soil, sleeper, etc) -which has to be best described.

Different numerical approaches are commonly used to model ballast. The choice amongst those strategies is usually done regarding the study scale of the track system ; switching generally between micro-and macro-mechanics approaches. Micromechanic approaches with Discrete Element Methods (DEM) are particularly well suitable to model the granular material behaviour of the ballast; whereas macromechanics approaches using Finite Element Methods (FEM) are preferred when the whole behaviour of the track system has to be described.

Those two approaches offer undoubtedly relevant simulation frames to perform analysis at different study scale. However, the fact remains at the same time that in both approaches drastic hypotheses are done to describe the system track regarding those frames. The influence of those hypotheses on the results have to be discussed. This is the issue of this work which aims to propose a new modelling srategy coupling FEM and DEM approaches to improve the description of ballast/soil interaction.

The analysis of ballast/soil interaction is traditionally performed using FEM. Ballast is therefore modelled as a simple layer of a continuous material. Using this modelling strategy, local behaviours due to the presence of a granular material in the track system -located in the layer of ballast (interaction network of grains, etc) or at the interface with the soil (phenomenon of stamping) -are not described.

As mentioned previously, DEM methods seem particularly well relevant to study such behaviours. Discrete models classically used in such studies take into account the polyhedral shape of the grains considered as rigid bodies. Such models make possible to perform rather realistic simulations of the ballast as a granular material. However, the choice of a poor bulk behaviour to model grains (rigid bodies) and the application of poor boundary conditions (rigid foundations) may imply a lack of deformability of the track system whose influence on its resulting behaviour has to be studied carefully by going further than introducing some relaxation in the model using a network of springs at the soil/ballast interface.

The originality of this work concerns the modelling approach based on the coupling between the standard DEM method and the FEM one. The long term purpose is to develop an accurate FEM/DEM railway track model making possible to provide more realistic responses of the system and better understand physic phenomena involved by studying the ballast behaviour itself and its interactions with the other track components (and with the soil especially).

One of the critical aspects we have to deal with when using DEM is the computational time which is strongly related to the number of bodies in the collection, the shape of the grains, the complexity of the bulk behaviours involved or coupling effects on contact problem in presence of large scale bodies. Simulations of railway track -even using poor grains modelling strategy -are from this point of view very expensive.

The strategies developed to succeed in performing accurate and efficient hybrid FEM/DEM simulations take place in the particular framework of the Non-Smooth Contact Dynamics (NSCD) method -a DEM method -proposed by Moreau and Jean. This framework is described in the first part of this paper. In the proposed hybrid model of the track system, grains are modelled as a collection of rigid bodies in interaction whereas layers of soil are described using finite elements. Three modelling options of the hybrid FEM/DEM model are discussed regarding the technical aspects of the coupling and the management of the contact approach more precisely, the physical behaviours which are possibly modelled with the proposed strategies and their efficiency.

The first approach makes possible to describe the whole local phenomena in the railway track both at the interaction soil / ballast and in the granular material (full micro-scale simulation). The second one is considered as a macro-description of the soil / ballast interaction to the extent that stamping effects of the grains in the soil are not modelled : due to an adapted contact approach developed in the open-source LMGC90 software used to perform the simulation, soil is modelled as a poorly discretised meshed finite element object whose characteristic length is much longer than the mean size of the ballast grains. If this approach is far less costly regarding CPU time, it does not allow to consider a full micro-scale analysis of the structure both inside the granular material and at the interface with the soil. That's why a prospective third approach is finally proposed to combine this new contact approach with the introduction of non linearity due to permanent plastic deformation of soil at contact locus thanks to a contact law miming the plastic behaviour and the local penetrations of the grains in the ground (stamping).

If the proposed FEM/DEM hybrid strategies aim to make possible the introduction of different components of the track system (soil, Under Sleeper Pad (USP), etc) as deforming bodies, the field of our study will be focused only on the influence of soil conditions on the ballast behaviour for different quasi-static load conditions. Further studies will have to be dedicated to the modelling of the other components as deformable bodies.

Framework of the Non-Smooth Contact Dynamics

DEM methods aims to model collection of bodies in interaction. Amongst DEM methods the distinguish is to be made between methods using explicit interaction laws defined with time derivative function -the " Smooth " methods -like Cundall method [START_REF] Cundall | A computer model for simulating progressive large scale movements of blocky rock systems[END_REF][START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF] and the " Non-Smooth " ones using on the contrary implicit laws (threshold laws, function with gap, etc) managing to take into account shocks between objects.

We will focus in this paper on one of those methods -the Non-Smooth Contact Dynamics (NSCD) method -proposed by Moreau and Jean [START_REF] Horkay | Unilateral contact and dry fruction in numerical simulation of deep drawing[END_REF][START_REF] Jean | Unilaterality and dry friction in the dynamic of rigid body collections[END_REF][START_REF] Jean | Frictional contact in collections of rigid and deformable bodies : numerical simulation of geomterials[END_REF][START_REF] Jean | The non-smooth contact dynamic method[END_REF][START_REF] Moreau | Unilateral contact and dry friction in finited freedom dynamics[END_REF][START_REF] Radjaï | Discrete numerical modeling of granular materials[END_REF] and imple-mented in the open-source LMGC90 software. Hereinafter are given the outline of the NSCD method which appears as a dedicated simulation framework to manage the modelling of both divided and continuous media as well as the strategy implemented in the LMGC90 software.

Mechanical and Simulation Framework

Multi-contact simulations aim to model the global behaviour of collections of solid bodies in interaction by considering both the motion of each body of the set of bodies and its interactions with its neighbours. The problem to be solved has therefore to be described using two sets of unknowns : global unknowns (or kinematic space unknowns) related to the bodies (center of inertia or mesh node displacement and velocity, resulting force and momentum, etc) and local unknowns (or contact space unknowns) related to interactions (gap between bodies, relative velocities, forces, etc). The idea is to switch from one description to the other one so as to be able to take into account both contributions of the interactions and the free motion of bodies in the global behaviour of the system.

Dynamics

For a standard formulation of the mechanic behaviour of deformable bodies, evolution equations of the dynamics of a set of bodies can be written in a generic form

Mq + K q + Cq = F ext (1)
M is the mass matrix, K and C are the rigidity and the viscosity matrix (with null contribution of rigid bodies), F ext are the external forces and q defines the associated degrees of freedom (the center of inertia or mesh node displacement here so that in the following q = V is a velocity). Since shocks, or other non-smooth events, can occur the acceleration and the contact forces must be viewed as differential measures. In the following, the integrated form over a time step will be used. It can be interpreted as a momentum balance over a time step h. Using a θ-scheme time integrator the Equation (1) can be written as :

(M+hθC+h 2 θ 2 K)(V i+1 -V i ) = h(1-θ)(F i +P i )+hθ(F i+1 +P i+1 )+I c, i+1 (2)
where F and P are the external and internal forces, and I c the impulsion due to contacts.

The behaviour of each body of the simulation is consequently described with a linear system of equations which can be written as :

M(V i+1 -V i ) = I f ree + I c, i+1 (3)
The resolution of the linear systems is performed using linear algebra library (Lapack, MUMPS, etc) or using local numerical methodologies implemented in the LMGC90 software. Equation (3) can be written in terms of velocities :

V i+1 = V i + M-1 I f ree V f ree + M-1 I c, i+1 (4)
where V f ree is the free (of interaction) velocity term.

At this step, the equations are described with global unknowns. The strategy to solve the problem consists in writing those equations in terms of local unknowns (at the contact scale) for which force/velocities relations are given by contact laws. In the following are given details surrounding the contact definition problem and the mechanism making possible to switch from one description scale to the other one. It has to be noticed also that this description makes possible to mix both rigid and deformable bodies in the same formulation.

Contact, Contactors and Contact-Body Mapping

If the system of bodies seems well defined according to the physical parameters (mass, Inertia matrix, etc), geometrical information allowing to describe the shape of the bodies are missing to totally describe the system. Those information are however essential to deal with the contact issue contained in the term I c of the Equation (4).

The management of the contact between bodies is based on the so-called slave/master approach or candidate/antagonist approach (Figure 1). According to this approach, a point C belonging to the well defined surface of contact with a candidate status of a body b may not go through an antagonist surface of a body a. For each situation of contact, a local frame (t, n, s) attached to the antagonist body is defined at the contact point A. Normal component of those frames n are directed towards the candidate body and correspond to the vector perpendicular to the tangential plane surface T defined by (t, s). In this local frame are described some local values (gap g, relative velocity V) inherent to the contact between the bodies in interaction.

Figure 1: Candidate/antagonist approach (2D and 3D).

Obviously, the management of the problem of contact needs to know the contact surface of each body in the simulation. Contactors play this role. They make possible to describe indeed the shape of the different bodies of the collection. Rigid contactors (disks/spheres, polygons/polyhedrons, cylinder, etc) are attached to the inertia center of the rigid particles and may play both roles of antagonist or candidate surface whereas contactors are defined for a given status -candidate or antagonist -through a set of surfacic nodes belonging to the meshed bodies when deformable. From this appears a natural mapping between local and global kinematic unknowns defined with a linear operator H(q) build on kinematic considerations (Equation ( 5)). Using duality considerations, a second linear operator T H(q) can be build (Equation ( 6)) so that the mappings between contact and body variables (i.e. relative velocities V and impulses I at contact and velocities V and resultant impulses I on body) is total (Equation ( 7)).

H(q) : (I α ) → I = α H α (q) I α (5)
Figure 2: Role of the linear operator H in the mapping of local and global velocities (contact between rigid and deformable bodies, and between rigid bodies).

T H(q) : V → (V α ) = α T H α (q) V (6) V ← Dynamics equation → I Body ↓ T H(q) ↑ H(q) V α ← Contact laws → I α Contact (7)
Thanks to these mappings and introducing the Delassus operator W = T H(q m ) M-1 H(q m ) for a configuration of contact q m , dynamics (Equation ( 4)) can be written in terms of contact unknowns:

V i+1 = V f ree + WI (8)
and be solved using contact laws defining the relation between relative velocities V and impulses I at contact for a couple of bodies in interaction.

It has to be noticed once more that the formulation proposed in Equation ( 8) is convenient for both rigid and deformable bodies. The framework proposed in the NSCD method seems hence suitable for performing hybrid FEM/DEM simulation.

Interaction Law

A classical velocity frictional contact law can be used to model accurately enough friction in granular material. The contact law described hereafter is currently used to model ballast grains interaction in railway track. Contact may be written as :

if g > 0, I n = 0 and any value for V n if g ≤ 0, V n ≥ 0, I n ≥ 0 and V n • I n = 0 (9)
and a Coulomb's friction law is used when g ≤ 0:

I t ≤ µ I n I t < µ • I n ⇒ V t = 0 I t = µ • I n ⇒ V t = -λ • It It ( 10
)
where λ > 0 and µ is the friction coefficient between the two bodies. This phenomenological law models dry friction between various materials but does not take into account local phenomena like local plastification or default abrasion; even if it can be enhanced by adding some local mechanisms like dilation or rolling friction. Many other interaction laws are available in LMGC90 making possible to take into account for example cohesion effect.

Resolution Strategies

Numerical Method

A contact α verifies the following set of equations :

     V α i+1 = V α f ree + α =β W αβ I β + W αα V α contact law (I α , V α ) = 0 (11)
It is possible to solve this set of equations if the impulses due to other contacts I β are known. This is the founding principle of the Non Linear Gauss-Seidel (NLGS) method used here. The contacts are solved one after another and iterations on all contacts are repeated until convergence. The algorithm is summed up below:

(0) Evaluating all the matrices W αβ           k = k + 1 (NLGS iteration)       α = α + 1 (Contact index) (a) Evaluating the right-hand side V α loc = V α f ree + β<α W αβ I k+1,β + β>α W αβ I k,β (b) Solving the local problem Convergence test for k = 0 . . . k max Table 1: NLGS algorithm.

Enhanced Strategies

Here are briefly described two numerical strategies -the multi-threading and the multi-processing approach -developped in the LMGC90 software to enhance computational performance.

An obvious method of parallelization is to use openMP directives around computation loops to obtain a shared memory multi-threaded code. This method is very easy to implement but needs tuning. It works exactly like the sequential method if the operations made inside the loops are independent. This method is used to parallelize the contact resolution loop. Even if the operations inside the loop are not strictly independent, the obtained results are good [START_REF] Renouf | Numerical investigations of fault propagation and forced-fold using a non smooth discrete element method[END_REF].

A domain decomposition method was also developed (Non Smooth Contact Domain Decomposition [START_REF] Visseq | Distributed Nonsmooth Contact Domain methodsDecomposition (NSCDD): algorithmic structure and scalability[END_REF]). One subdivides the domain and solves the contact dynamic problem in each sub-domain with NLGS algorithm; gluing conditions ensure velocity continuity of interface bodies (when converged), like in FETI method [START_REF] Farhat | A transient FETI methodology for large-scale parallel implicit computations in structural mechanics[END_REF]. This strategy can be used only for a limited set of rigid contactors (sphere, polyhedrons).

Those strategy undoubtedly enhance the computational performance of the LMGC90 software. An effort has to be made to extend the NSCDD method to the whole set of contactors. This method can not be used yet to perform simulation of hybrid FEM/DEM model. Regarding multi-threading strategy usable with hybrid models, it have not been used to perform simulations presented hereafter. Further studies could explore the impact of this approach on computational times.

Description of a Simplified Track System and Modelling Strategies

Simulation of hybrid models with different modelling strategies are presented in Sections 4, 5, and 6. Before focusing on simulation results, this section aims to clearly fix and justify the modelling hypotheses used in the following work (geometrical description of the samples, bulk behaviours of the different components of the railway track, loading scenarios, etc).

Description of the railway track system

Basically track system may be defined as a succession of layers associated to the different components (soils, ballast, USP, selleper) of the track which forms embankments whose profiles may slightly differ from one to another track. A simplified description of the track is proposed here. It consists mainly in considering a multi-layers square portion of track as shown in Figure 4 with periodic conditions so as to model the track as a continuous media. A polyhedral sleeper (h sleeper x L sleeper x L sleeper ) is anchored in a ballast bed (h ballast , lateral anchorage L anchorage ) which sits on top of layer(s) of soil(s) (h soil(s) ). USP is not modelled in this work. Dimensions may change for the different presented simulations. Precisions are given when needed. Simulations are performed in the following with a single layer of soil.

Modelling Strategy

Ballast is modelled as a coarse granular material using the DEM approach. Grains are considered as rigid bodies (ρ grains = 2700 kg/m 3 ) with frictional behaviour (Section 2.1.3) to model the contact between grains (µ grains = 0, 80).

The soil is modelled as a deformable body using the FEM approach. A basic elastic linear behaviour is chosen to describe the response of the layer of soil whose mechanical characteristics are summed up hereafter : ρ soil = 2000 kg/m 3 , Young modulus E soil to be set, ν soil = 0, 3. This choice has been made to reduce as far as possible computational time; the use of non linear bulk behaviour models implying over calculation costs (matrix of rigidity K and viscosity C evolving and having to be rebuild during all the computation at each time step k (K k and C k )). More complex bulk behaviours (plasticity particularly) should be tested to model the layer of soil. The contact between soil and ballast is modelled with friction too (µ grains/soil = 0, 80).

Considering the discretisation of the meshed layer of soil, it is highly related to the choice of the contact strategy to be used in the simulation. The choice of the type of contactors -antagonist or candidate -affected to both grains and the upper face of the soil column to describe the contact between the ballast and the ground may impose constraints on the mesh and on the definition of the characteristic length of the surfacic mesh on which grains are in contact particularly. Two strategies have to be considered depending on whether the upper surface of the column of soil plays the role of candidate or antagonist surface (Figure 5). On the one hand, the upper surface of the layer of soil may have a candidate status. Hence the discretisation of the contact locus depends on the mesh discretisation. This is due to the fact that the possible locus of contact points are in this case either confounding with the nodes of the mesh either located at the point projection of the Gauss point of each surfacic element of the mesh. Therefore the mesh discretisation has to be adapted to the size of the grains of ballast so that the contact between the ballast and the soil is accurate (no interpenetrations). This is the strategy used in Section 4 which makes possible to model local phenomenologies at the ballast/soil interface like stamping. On the other hand, the free surface of soil in contact with grains may be described as an antagonist surface. In this case, the discretisation of the soil has no influence anymore in the management of the contact between the ballast and the soil to the extent that contact locus are directly carried by the contactors defining the shape of the grains. This strategy is used in the simulations presented in Section 5. If this strategy is all the more efficient that the mesh is roughly discretised, it makes possible to catch a macroscopic soil/ballast interaction behaviour only. A very last strategy presented in Section 6 aims to counterbalance this aspect by adding a specific contact law in the model so as to manage some plastification at the contact scale. Strategies are summed up in Figure 6. The limitations of those modelling strategies are discussed in detail hereafter regarding the phenomenologies they may model and their needs in computer ressources to perform simulations.

Loadings and Boundary Limits Conditions

Nodes belonging to the lower surface of the column of soil are fixed (null velocity applied in the 3 space directions) whereas the upper contact surface is free of motion. Lateral surfaces of the column of soil are also not constrained. No particular needs in numerical strategy -like absorbing boundary conditions for example -to deal with wave reflection have been noticed. This is mainly due to the type of loads applied (quasi-static) during our simulation. Some dissipation is introduced in the model throughout the θ-scheme (θ slightly higher than 0.5) and Rayleigh damping (α and β equal to 7%).

Regarding loading conditions, two main types of loads have to be distinguished. Compressive vertical actions are applied on the sleeper using punctual forces F v . The maximal magnitude of the force is defined so that the stress under sleeper σ v reaches 100 or 150 kPa. This compressive action is progressively applied on the center of inertia of the sleeper. Lateral actions may be applied on the sleeper too. Velocity in the x direction V x is applied to the center of inertia of the sleeper so as to obtained a cumulative displacement of the sleeper of 5 mm at the end of the simulation. Of course, those different loadings may be combined as shown in Figure 7 (a) which describe loading actions for a resistance test to the lateral displacement under vertical loading. We may also apply a succession of loading and unloading actions to the sleeper to perform endurance tests in a range of given frequencies (Figure 7

Overview of the Interest of the FEM/DEM Coupling

Even if this paper is mainly devoted to the strategies used for hybrid modelling some general results are given in the following to show the interest of the coupling in the simulation of railway track. If both macroscopic and microscopic analyses may be performed to measure the role of the deforming foundation on the response of the railway track, we will focus here on the description of the global behaviour of the track system in presence or not of deforming foundation throughout very simple and well known indicators. One may describe for example the evolution of the settlement of the sleeper over the time or the force-displacement response of the sleeper as in Figure 8 for which sleeper of a 60x60 cm 2 sample is submitted to the application of loading and unloading cycles of a vertical load F v as described in the endurance test Figure 7 In this case, one may measure the role played by the soil in the global response of the track system. Energy is stored in the deforming layer of soil and then returned back to the system during unloading phasis (relaxation). This is not the case when the foundation is modelled as a rigid body. Global response seems hence to be more conservative. Of course those very first FEM/DEM responses might be considered with a lot of caution : complete qualitative parametric studies have to be performed so as to definitely asses the pertinence of using hybrid models in railway simulations. Inhere we will stay focus nevertheless on the adopted strategy to model systems with the coupling approach.

A Full Microscale Modelling Option to Describe Local Phenomena

The very first strategy which is proposed in this paper consists in a full microscopic approach of the modelling of the railway track system both inside the granular material and at the interface between soil and grains. If the scale of the analysis of the layer of ballast is frozen due to the use of DEM to model ballast as a coarse material, the way to deal with the discretisation of the layer of soil makes possible to model local phenomenologies at the boudary of the FEM and the DEM models ; the choice of the discretisation being related itself to the contact strategy defined to manage the soil/ballast interaction as mentioned in Section 3.2. In this approach, the upper surfacic mesh of the column of soil carries candidate contactors so that the characteristic length of the mesh at the boundary has to be strictly related to the size of the grains in contact to manage the contact detection. Regarding our study case, it implies to perform the mesh of the layer of soil imposing a maximal size of the element of 1 cm, which nearly corresponds to the minimal radius of the sphere in which the smallest grains of ballast can be encapsulated. This condition on the meshing process has to be obeyed because otherwise contact detection failure may surely happen with grains going through the soil of layer as if there were no ground at all (Figure 9).

Meshing Strategies of the Column of Soil

If the use of the finest discretisation may prevent problems of detection, it takles other issues espacially concerning the management of the computing needs and the control of computational costs. Those issues are all the more prevailing than the size of the railway track to be modelled is important. Of course, the size of the corresponding linear system to be solved may strongly expand with the dimensional characteristics of the sample (increase of the number of Degrees of Fredom DOFs of the system) as well as the number of coupled contacts (H(q) and T H(q)operators). Different meshing strategies have been considered in this work so as to reduce as much as possible the needs imposed by the management strategy of the problem of contact. Anisotropic meshes in the z diection are currenly used so as to build the most economical system in a computational point of view. The size of the surfacic elements of the upper face of the column of soil respects the criterion imposed to well manage detection of contact and insure the non interpenetration of grains in the soil whereas a rough discretisation is used to describe the bottom of the foundation (20 cm versus 1 cm). This meshing strategy makes possible to considerably reduce the number of DOFs and the number of volumic elements in the system whilst insuring a suitable contact detection. Figure 10 shows the impact of this strategy in comparison to different meshes of a column of soil with a varying length side of its square section (L sample between 30 up to 200 cm) and a constant height (h soil = 40 cm). Meshes are performed with the free Gmsh software [START_REF] Geuzaine | Gmsh : a three-dimensional finite element mesh generator with built-in pre-and post-processing facilities[END_REF]. and meshing strategies (homogeneous mesh with characteristic length equal to 1 or 2 cm, anisotropic mesh in the z direction).

Description of Local Phenomena like Stamping

Once the difficulties encountered to set the microscale model managed, results of simulation make possible to describe very local phenomena. Stamping effetcs due to the penetration of grains of ballast in the layer of soil are well described as shown in Figure 11 (a) which represents an hybrid FEM/DEM sample (30x30 cm 2 ) under a vertical compressive load applied throughout a plane. The punctual load is increased up to reach 5 T which is not representative regarding railway track loading actions. Such a compressive load has been applied in the sample so as to highlight as much as possible the local phenomenologies which can be catched using this modelling approach.

In the vertical cross-section of the sample Figure 11 (b), blobs with a high level of stress may be located at the soil/ballast interface. Those local effects affect a thin transition layer of the ground which returns to an homogeneous stress level in the layers below. Deformations of the foundation due to the grains is shown Figure 11 (c). Grains marks are visible and form a speckle pattern.

This type of simulation may be considered only for small samples like the one presented herebefore. Even if results give a good description of the behaviour at the soil/ballast interface, computational costs are prohibitive. The load of the 30x30 cm 2 sample has needed 17 hours of computing with an anisotropic mesh of the foundation (characteristic length of the mesh at the surface of contact grains/soil equal to 2 cm) while the model is made of only 254 rigid grains of ballast. Computing delay exploses with a finer discretisation, reaching 170 hours of CPU for a characteristic length of the mesh of 1 cm. The needs in memory to deal with the simulation and to load the system to be solved are another blocking point of the microscale strategy. The use of such a strategy may not be used for realistic samples due to the huge number of DOFs envolved to model the soil and or coupled contacts at the interface soil/ballast. Modelling difficulties appear for simplified sample of the railway track even for modest sizes (45x45 cm 2 square cross-section).

All those difficulties encountered regarding the performance of the microscale models lead to prefer other strategies of simulation to deal with soil/ballast interaction issue as the following ones presenting in Sections 5 and 6.

Toward a Macroscale Model of the Soil/Ballast Interaction

Managing realistic simulations of the track system has been the very purpose of the strategy presented in this section. First approach has shown the difficulty to both deal with the technical management of the problem of contact -and the detection between the grains and the soil especially -and the modelling of local behaviours for realistic systems of track whose dimensions are much more important than the one dealt with in the simplified study cases. A typical realistic slice of track which could be studied for example is made of a 4,5 m length, 1,0 m large and 50 cm heigth layer of ballast composed of nearly 30000 grains (Figure 12). Regarding the geometrical size of such a system, a new approach of the contact detection has been considered so as to be able to deal with such simulations in terms of efficiency : insuring a good contact detection for a controled CPU time and memory usage has been the challenges to meet with the management of the macroscale approach. Herein the interaction of the soil and the ballast is poorly described. This is the reason why the strategy proposed in this section is the so-called macroscale approach. The idea is to perform the mesh discretisation independently of size considerations -characteristic length of the meshing being to be set regarding the size of the ballast grains -imposed by the strategy of contact contact used for the detection. The opposite contact detection strategy is used in this approach which means in contactors terms that the upper surfacic surface of the column of soil has an antagonist status instead of a candidate one. From this point of view, grains which are now considered throughout candidate contactors may not go through the meshed column of soil independantly of the size of the surfacic elements at the contact interface. This approach makes possible to deal with homogeneous hexaedric mesh with a limited number of elements. Hence a limited number of DOFs is used to describe the soil behaviour.

Performance and Reliability

If gains in terms of performance are indisputable for the roughest discretisation (Figure 14), it remains true that the model may not permit to describe local phenomenologies (Figure 13 (b)). Stamping for example may not be modelled using this approach except using a very fine discretisation of the column of soil ; which appears as a non natural usage of the proposed strategy to be undoubtedly proscribed. Nevertheless, in this condition, local phenomena may be described in the same way as in section 4 as shown in the comparative Figure 15. In the same time it remains that the strategy is no more efficient regarding the reductions of CPU time and the memory usage during the simulation. Simulation may be performed perhaps for an higher computational cost due to the nature of the volumic elements composing the mesh. The proposed strategy may or may be not consider as relevant to model realistic systems depending of the point of view. If the success of the strategy has to be considered throughout efficiency criteria, the challenge is clearly met : the macroscale approach makes possible to manage heavy samples. But it is all the contrary if the finality of the simulation lies in modelling local phenomena like stamping presented in Section 4. In this case choices to be made to keep a saving time strategy are incompatible with details to be catched in those types of behaviours. 6 From Micro to Macroscale Modelling : On the Use of a Specific Contact Law to counterbalance : the "Plastic" Law

The very last modelling approach aims to go further in the simulation by both dealing with realistic track systems and providing a local response of the model at the interface between soil and grains. To our minds the strategy developped in Section 5 is the strategy to be used in the case of hybrid FEM/DEM modelling of railway systems so as to keep performance gains. But tkis strategy has to be adapted to take into account local effects of the ballast on the deforming foundation. This is the role to be played by a new interaction law which aims to couterbalance the lack of micro mechanisms at the contact scale between grains and soil by introducting some yielding effects and a controlled interpenetrations of the grains in the column of soil. The so-called "plastic" law is detailled in the following as well as some results of simulation in which the law is used to show its influence. This work is very prospective and may have to be enhanced before beding used for simulations taking place in an engineering study for example.

The "Plastic" Interaction Law

The plastic law aims to model friction but differs from the philosophy of the unilarity condition of Signori presented in Section 2.1.3 to the extent that it manages interpenetrations of objects. Nevertheless the law is very closed in its formulation either considering the normal or the tangential behaviour. Friction behaviours is slightly modified in comparison with the classical Coulomb law to take into account dilantancy phenomenon due to the interpenetrations.

Normal Response :

If the unilarity condition traditionaly describe the fact that a body may not go through another one as mentioned in Section 2.1.3, the plastic law allows some interpenetrations g of the objetcs in contact. In the case of the plastic law, the formulation of the normal component is only differing in the condition on the gap g between objects, which has not to be positive or nul anymore as in the classical writing (Equation 9) but has to be greater or equal to a maximal autorised gap ḡ. Thus the condition of normal contact may be writen as : if g > ḡ, I n = 0 and any value for

V n if g ≤ ḡ, V n ≥ 0, I n ≥ 0 and V n • I n = 0 (12)
A graphical representation of the unilarity condition taking into account autorised interpenetration ḡ is over ploted in Figure 16 (a The interpenetration is managed thanks to an evolution law (Figure 16 (b)) which makes possible to define the value of the autorised interpenetration ḡ during the computation. If the normal reaction I n is greater than a critical value I crit n set in the parameters of the law, then the interpenetration of bodies is allowed and the ḡ parameter is incremented of dḡ (an evolution parameter (dḡ > 0) set in the definition of the law too). In this way, the interpenetration between two objects goes incrementaly from a nul interpenetration to a maximal value ḡlim equals to 20 % of the radius of the encapsuling sphere of the candidate body.

A sort of hardening phenonmenon is also introduced : the value of the critical reaction from which depends the stamping mechanism being defined regarding the value of the parameter ḡ : I Friction gives a relation between the tangential reaction I t and the sliding velocity V t . Generally, contact law between bodies model dry friction, the most widespread being the Coulomb law as described in Equation 10.

In the case of the plastic law, the friction coefficient µ(ḡ) change when interpenetration occurs going from µ s to µ p using an increment dµ. µ s and µ p are the static and plastic friction coefficients defined in the parameters of the plastic law as well as the second evolution parameter dµ. The expression of µ(ḡ) is given by : µ(ḡ) = min(µ s + N × dµ, µ p ) with N positive or integer. The increase of the friction coefficient aims to take into account dilatancy effects implied by stamping.

Influence on the Relaxation of Samples

Using this contact law and the macro strategy in the contact detection is a way to introduce relaxation in the simulation. Full microscale models at the interaction scale makes possible to naturally relax the granular system thanks to the local rolling of grains and the local stamping at the locus of contacts and even more to catch the effects of the stamping in the response of the soil layer (blobs/peak stress). This is not exactly the case with the strategy presented herein to the extent that the specific contact law used may only introduce some relaxation (by managing interpenetrations at the contact between grains and foundation) without acting on the mchanical behaviour of the deforming layer of soil. The response of the soil depending on the mesh discretisation remains still very macroscopic. To be convinced of this, we may consider the response of a simple FEM/DEM sample made of 4 cubic rigid bodies with different sizes (c 1 = 25 cm, c 2 = 5 cm, c 3 = 12 cm, c 4 = 30 cm) under self weight settled on a deforming foundation modelled using the poorest discretisation ever (a single hexaedric element) as depicted in Figure 17. Contact between the blocks and the foundation is modelled with the plastic law without hardening. The same parameters to define the law are used so that the settlement is only depending on the size of the blocks. At the end of the simulation, blocks have reached the maximal autorised interpenetration -id est 20% of c i , i in [START_REF] Cundall | A computer model for simulating progressive large scale movements of blocky rock systems[END_REF][START_REF] Geuzaine | Gmsh : a three-dimensional finite element mesh generator with built-in pre-and post-processing facilities[END_REF] -, which shows that "geometrical" stamping is currently modelled whereas the response of the soil remains macroscopic.

Even if the local influence of the stamping on the response of the deforming column of soil is not described with this stratgy, effects of the geometrical relaxation of the layer of ballast are clearly perceptible on the response of the sample either at the FEM/DEM interface or in the granular material as shown with the 60x60 cm 2 simplified sample whose sleeper is loaded over time (F v ). Results for two discretisations of the layer of soil (5x5x2 elements in the x, y and z directions or 25x25x10) are given Figure 18 at the soil/ballast interface (a) or in the granular material itself (b). Plastic law carries out a more homogeneous contact distribution inside the granular material or at the boundary with the soil. Consequently, strong interactions networks are broken due to local grains rearrangements that this strategy manage to model.

Here are given qualitative results. This modelling strategy has to be tested more carrefully. Several issues have to be tackled. The very first one to be considered is the setting of the parameters of the law. Parameters used in the simulations presented in this Section have been chosen with purely numerical purpose. A physical interpretation has to be given more especially. This tackle another issue : the formulation of the proposed law and the definition of its evolution scheme, which seem very poor. A modelling effort seems to be necessary to improve the contact law before using this law with the macrosclae modelling strategy of the soil/ballast interface to model railway systems efficiently and in a more accurate way.

In this paper three hybrid FEM/DEM modelling strategies of the soil/ballast interface are described using the NSCD framework. Those modelling strategies make possible to take into account the behaviour of the soil/ballast interface in more or less details. Very fine descriptions of phenomenologies may be done using the full microscopic approach. The mechanical behaviour of the soil and of the soil/ballast interface are then well described when stamping occurs. The main drawback of this type of modelling strategy remains in the huge needs in computing ressources accompanying which go along with. More economical approaches to be used are presented as well taking into account with more or less accuracy the local phenonmenologies involved in the soil/ballast interface. The second method which focuses only on computing performance does not take into account those local behaviours at all whereas the third one aims to introduce those phenomenologies using a specific interaction law -whose management has to be improved -to well describe the grain/ballast contact. The three approaches are technically working and may be used for the modelling of railway systems. The choice of a given approach rather than another one is done at the discretion of the modeller depending on the phenomenologies to be modelled in the simulation and/or the computing ressources.

If the feasability of the hybrid FEM/DEM coupling in the railway framework has been shown, further work has to be done nevertheless to tackle some modelling issues which have not been discussed in this paper. It mainly concerns the FEM description of the colmumn of soil. A linear elastic behaviour is used on a poorly discretised column of soil to keep saving CPU time. Of course these choices have to be discussed by taking into account first of all mesh convergence of the response of the soil in our calculations and by using more realistic non linear bulk behaviour -introducing yielding phenomena for example -to describe the column of soil. Moreover no particular numerical treatments have been carried out to manage wave propagation/reflection in the deforming layer. This can be explained throughout the loads considered in our simplified study cases. This is a key point which has to be studied so as to be able to use hybrid FEM/DEM with dynamic sollicitations.

Considering more realistic sollicitation in the simulation has to be done while taking into account more realistic description of the samples too. This implies to be able to deal with samples of important size. This is another key point we have to meet with to measure the influence of the FEM/DEM coupling on characteristics responses of railway systems. The qualitative study work which has been done herein has to be completed with quantitative research. This work seems inseparable with the use of numerical strategies making possible to enhance computing performance. Multiprocessing methode NSCDD has to be extended in this way to be able to deal with hybrid FEM/DEM sample together with or not multithreading. In this manner, hybrid models will be defininitely adopted in the modelling of railway track systems or not. In such a case, the coupling upon one of the suitable strategy may be used to optimize engineering common restoration procedures like tamping process for example.

Figure 3 :

 3 Figure 3: Signorini-Coulomb law: friction with unilaterality condition.

Figure 4 :

 4 Figure 4: Simplified geometry of the track system.

Figure 5 :

 5 Figure 5: Contact strategies for the ballast/soil interaction : (a) upper surface of the soil as candidate surface / grains of ballast as antagonist ones, (b) inversion of the status of contactors.

Figure 6 :

 6 Figure 6: Modelling strategies and overview of the described phenomenologies.

Figure 7 :

 7 Figure 7: Loading actions for a mixed vertical and horizontal sollicitations test (a) or for endurance tests with variation of the sollicitation frequency (b).

  (b).

Figure 8 :

 8 Figure 8: Influence of deforming foundation (in black) on the macroscopic response of the railway track (comparison to the full rigid case in red) : (a) time evolution of the settlement and (b) Force-Displacement response of the sleeper submitted to 5 compressive loading cycles (frequency f=5 Hz (a)) and respectively 1 cycle (frequency f=20 Hz (a)).

Figure 9 :

 9 Figure 9: Mesh discretisation and problem of contact detection for a characteristic length of the mesh equal to 2 cm : (a) overview of the layer of ballast relying of the deforming foundation and (b) detail of grains going through the meshed column of soil.

Figure 10 :

 10 Figure 10: Mesh characteristics : (a) number of DOFs and volumic elements for different length side of a square cross-section column of soil (constant height), (b) overview of the mesh for different sizes of the sample (30x30, 45x45 and 75x75 cm 2)and meshing strategies (homogeneous mesh with characteristic length equal to 1 or 2 cm, anisotropic mesh in the z direction).

  Figure 10: Mesh characteristics : (a) number of DOFs and volumic elements for different length side of a square cross-section column of soil (constant height), (b) overview of the mesh for different sizes of the sample (30x30, 45x45 and 75x75 cm 2)and meshing strategies (homogeneous mesh with characteristic length equal to 1 or 2 cm, anisotropic mesh in the z direction).

Figure 11 :

 11 Figure 11: Microscopic approach for the modelling of stamping effects of ballast grains on the upper surface of the layer of soil : (a) 3D sample, (b) vertical crosssection with identification of local blobs over stressed and (c) deforming horizontal cross-section of the interface soil/ballast (shape marks of grains).

5. 1 AFigure 12 :

 112 Figure 12: A strategy to deal with the modelling realistic railway track : (a) railway track and (b) representative hybrid FEM/DEM model of a slice of the track.

Figure 13 :

 13 Figure 13: Influence of the inversion of the contact detection strategy on the described phenomenologies : (a) simplified hybrid FEM/DEM sample with the macroscale approach under self weight with developped interaction network and (b) normal stress on the deformed column of soil with dependency to the mesh.

Figure 14 :

 14 Figure 14: Influence of the discretisation of the column of soil on CPU times for the simulation of the column under sef weight (a) and the simulation of the whole track system (soil + column) under self weight and under a vertical load F v applied to the inertia center of the sleeper (b).

Figure 15 :

 15 Figure 15: Comparison of horizontal cross-section for micro-(a) and macro-(b) approaches for similar characteristic length of the meshes.

Figure 16 :

 16 Figure 16: Introduction of "plastic" phenomena at the interaction scale using a specific interaction law managing a controled interpenetration : (a) unilarity condition with the maximal autorised interpenetration ḡ, (b) evolution function managing the penetration of the candidate in the antagonist contactor (c) tangential behaviour modelling pseudo dilatancy.

  crit n (ḡ) = I crit,0 n + ḡ × k where k corresponds to the hardening coefficient and I crit,0 n the critical value of the local reaction. Those two parameters are needed in the definition of the law. Tangential Response :

Figure 17 :Figure 18 :

 1718 Figure 17: Settlement of rigid blocks in a deforming foundation : macroscopic response of the foundation (a) and local settlement of the blocks : vertical interpenetration of blocks monitored over time (b).