Member, IEEE Claudio Weidmann
email: claudio.weidmann@ieee.org

Senior Member, IEEE Michel Kieffer
email: .kieffer@lss.supelec.fr

Evaluation of the Distance Spectrum of Variable-Length Finite-State Codes

Keywords: Variable length codes, finite state machines, source coding, channel coding, communication system performance. I

The class of variable-length finite-state joint sourcechannel codes is defined and a polynomial complexity algorithm for the evaluation of their distance spectrum presented. Issues in truncating the spectrum to a finite number of (possibly approximate) terms are discussed and illustrated by experimental results.

that the rate (output bits per input symbol) is not assumed constant per transition in the state diagram. Thus it is not possible to enumerate the information (source) spectrum together with the code spectrum; only the latter will be computed.

The contributions of this letter are as follows. Section II introduces variable-length finite-state codes and computes stationary state probabilities adjusted for variable-length transitions, which are needed to average the spectrum, improving upon the approximate approach in [START_REF] Ben-Jamaa | Analytical tools for optimizing the error correction performance of arithmetic codes[END_REF]. Section III presents the union bound and the spectrum, Section IV gives a recursive spectrum evaluation algorithm, while Section V presents a direct algorithm using numerical matrix inversion. Section VI discusses the issues involved in bounding the error when approximating the spectrum by a finite number of terms, which are themselves approximations. Finally, Section VII provides experimental results for an example JSC-AC.

II. VARIABLE-LENGTH FINITE STATE CODES

A binary finite-state encoder (FSE) may be represented as a directed graph Γ defined by a set S of states (vertices) and a set T of transitions (directed edges) between states, where each transition u ∈ T is labeled with a vector of input symbols I(u) and a vector of output bits O(u), whose lengths are ℓ(I(u)) and ℓ(O(u)), respectively. In the variable-length (VL) FSEs considered here, both input and output labels may be of variable length, generalizing fixed-length (FL) in -FL out FSEs [START_REF] Pollara | Finite-state codes[END_REF] and FL in -VL out FSEs [START_REF] Piret | Comma free error correcting codes of variable length, generated by finite-state encoders[END_REF]. The starting point is a FSE with nonempty input and output labels, e.g., a VLC for a K-ary source will be represented by the root state and K VL-output transitions, while a JSC-AC will be represented by a reduced finite-state automaton [START_REF] Ben-Jamaa | Analytical tools for optimizing the error correction performance of arithmetic codes[END_REF] with VL input and VL output labels, all nonempty. Since these are JSC encoders, each transition will also have an associated transition probability P (u), which may be computed using the source model. The sum of the outgoing transition probabilities of a state must equal one. In the following, a discrete memoryless source (DMS) is assumed, for which the transition probability is simply the product of the probabilities of the input label symbols. For a FSE to be a proper source encoder, for every state, the input labels of the outgoing transitions have to form a complete prefix set, which also implies that their transition probabilities sum to one.

Given an initial state s 0 , the encoding of all possible (semi-) infinite input sequences can be displayed by a trellis with states aligned at equal output bit length connected by VL transitions. The VL finite-state code (FSC) C(Γ, s 0 ) is the 0000-0000/00$00.00 c 2010 IEEE set of sequences of concatenated output labels corresponding to all paths through the trellis. Let σ(u) be the originating state of a transition u ∈ T and τ (u) its target state. A path

u = (u 1 • u 2 • . . . • u k ) ∈ T k on the trellis is a concatenation (denoted by •) of transitions that satisfy σ(u i+1 ) = τ (u i ) for 1 ≤ i < k (this corresponds to a walk of length k on Γ). The probability of a path is P (u) = k i=1 P (u i ).
For the purpose of computing the code spectrum, it is advantageous to consider a bit-clock trellis [START_REF] Ben-Jamaa | Analytical tools for optimizing the error correction performance of arithmetic codes[END_REF], where each transition is labeled with exactly one output bit. The corresponding bit-clock FSE Γ b (S b , T b ) is obtained by inserting additional deterministic transitions and states in order to split up output labels longer than one bit. A transition t with such a label is replaced by a chain of l = ℓ(O(t)) transitions t (1) , . . . , t (l) (and l -1 deterministic intermediate states), each labeled with the corresponding output bit, where t (1) "inherits" the input label of t and the corresponding probability, while t (2) , . . . , t (l) have empty input labels and probability one. Define M = |S|, then the state set S b ⊃ S of Γ b will be of size

M b = |S b | = M -|T | + t∈T ℓ(O(t)).
The following definitions will be used in the sequel. The set of transitions from x to y is

T xy = {u ∈ T b : σ(u) = x, τ (u) = y}, while T * (y, z) = {(u, v) ∈ T 2 b : σ(u) = σ(v), τ (u) = y, τ (v) =
z} is the set of all pairs of transitions ending in y and z, respectively, and not starting from a common state. The length of a path will exclusively denote its output length in bits, which equals the number of transitions on the bit-clock trellis. The Hamming distance between the outputs of a pair of n-bit paths (u, v) is denoted by d H (u, v). The set P n (x, y, z) contains all pairs of paths of length n bits starting from state x and ending in y and z, respectively, without merging in between, i.e., pairs of paths that have no state in common at bit indices 2, . . . , n-1. The set of all pairs of paths diverging in the starting state x and converging for the first time n bits later is

C n (x) = y∈S P n (x, y, y).
It will be assumed that the encoder graph is irreducible,1 i.e., that any state can be reached from any other in a finite number of transitions, and that it is aperiodic,2 i.e., the state recurrence times are not multiples of an integer period m > 1. These assumptions imply that the FSE forms an ergodic Markov chain, which has a unique stationary state distribution. Thus the asymptotic steady-state performance for long (unterminated) code sequences will be independent of the initial state. It will be weighted with probabilities p(x), x ∈ S, that the current code bit belongs to a transition leaving x, since paths can only diverge in states x ∈ S. Hence the probabilities p(x) will depend on the lengths of the transition labels in Γ.

Proposition 1: The stationary probability that the current code bit belongs to a transition leaving x ∈ S is

p(x) = p * x 1 + y∈S l≥2 (l -1)p l xy x∈S p * x 1 + y∈S l≥2 (l -1)p l xy , (1)
where p * x is the stationary state probability on the FSE graph Γ,

p * x = y∈S p * y   l≥1 p l yx   (x ∈ S),
and p l xy = u∈Txy:ℓ(O(u))=l P (u) is the probability of reaching y ∈ S from x ∈ S by emitting an l-bit output label.

Proof: By its definition, p(x) is the sum of the stationary probabilities (on the bit-clock graph Γ b ) of x and of all intermediate states x ′ ∈ S b \ S corresponding to transitions leaving x. Let q x and q l,k xy , respectively, denote these probabilities, where the latter stands for the k-th intermediate state of an l-bit transition from x to y (parallel transitions of the same length are merged and their probabilities added). The equations for the stationary distribution on Γ b are thus

q x = y∈S q y p 1 yx + l≥2 y∈S q l,l-1 yx (x ∈ S) (2)
q l,1 yx = q y p l xy (x, y ∈ S) (3)
q l,k yx = q l,k-1 yx (x, y ∈ S, 2 ≤ k < l). (4)
Combining ( 3) and (4) yields q l,l-1 yx = q y p l xy , which inserted into (2) leads to (1) after normalization.

III. UNION BOUND AND DISTANCE SPECTRUM

Consider a ML decoder applied to a FSC transmitted over a binary-input memoryless channel and let P e be the event error probability at any position in the code sequence. The following union upper bound holds [START_REF] Viterbi | Principles of Digital Communication and Coding[END_REF]Ch. 4.4]:

P e ≤ ∞ d=dfree A d P d , (5)
where

d free = min n≥1,x∈S min (u,v)∈Cn(x) d H (u, v) (6)
is the free distance of the code,

A d = x∈S p(x) ∞ n=d (u,v)∈Cn(x) dH(u,v)=d P (u) (7)
is the spectral coefficient that counts the average number of paths at Hamming distance d from a given path, where p(x) is the stationary probability of diverging in state x, and P d is the probability that the decoder selects an erroneous path at distance d instead of the correct path. For BPSK signaling with energy E b per channel bit over an AWGN channel with zero-mean Gaussian noise of variance N 0 /2 this term is [1, Ch. 2.9]

P d = 1 2 erfc d E b N0 ≤ exp -d E b N0 (8)
An important fact is that the union bound ( 5) is also a bound on the average bit error rate in the code domain, because an event error causing d code bit errors has length at least d bits.

IV. RECURSIVE SPECTRUM EVALUATION

Let A x d,n (y, z) be the average number of pairs of n-bit paths at Hamming distance d, starting in state x ∈ S and ending in states y ∈ S b and z ∈ S b , respectively, without merging in between,

A x d,n (y, z) = (u,v)∈Pn(x,y,z) dH(u,v)=d P (u). (9)
When z = y, paths merge for the first time at length n. Using (9), one may rewrite [START_REF] Jaspar | Design and performance analysis of joint source-channel turbo schemes with variable length codes[END_REF] as

A d = x∈S p(x) ∞ n=d y∈S A x d,n (y, y). ( 10
)
Proposition 2: For any x ∈ S, y, z ∈ S b , A x d,n+1 (y, z) may be evaluated recursively as

A x d,n+1 (y, z) = (u,v)∈T * (y,z) dH(u,v)=0 P (u)A x d,n (σ(u), σ(v)) + (u,v)∈T * (y,z) dH(u,v)=1 P (u)A x d-1,n (σ(u), σ(v)), (11)
for d ≥ 1 and n ≥ d. For d = 0 and n ≥ 1, only the first sum in (11) is computed (i.e., A x -1,n = 0 by definition). The recursion is initialized with

A x d,1 (y, z) = (u,v)∈Txy×Txz dH(u,v)=d P (u) (d = 0, 1), (12)
Proof: The initialization (12) is a simple reformulation of (9) for d = 0, 1 and paths of length n = 1 bit. To prove (11), one may write [START_REF] Piret | Comma free error correcting codes of variable length, generated by finite-state encoders[END_REF] for pairs of paths of length n + 1,

A x d,n+1 (y, z) = (u,v)∈Pn+1(x,y,z) dH(u,v)=d P (u). (13)
Consider any pair of paths (u, v) ∈ P n+1 (x, y, z). Since the two paths have not merged in between, there exist a pair of states (y ′ , z ′ ), y ′ = z ′ , a pair of paths (u ′ , v ′ ) ∈ P n (x, y ′ , z ′ ), and a pair of transitions (u, v) ∈ T * (y, z)

such that u = u ′ •u and v = v ′ •v. Two cases have to be considered. First, one may have d H (u, v) = 1, in which case d H (u ′ , v ′ ) = d -1. Second, one may have d H (u, v) = 0 and d H (u ′ , v ′ ) = d.
Hence (13) may be rewritten as

A x d,n+1 (y,z) = (u,v)∈T * (y,z) dH(u,v)=0 (u ′ ,v ′ )∈Pn(x,σ(u),σ(v)) dH(u ′ ,v ′ )=d P (u ′ )P (u) + (u,v)∈T * (y,z) dH(u,v)=1 (u ′ ,v ′ )∈Pn(x,σ(u),σ(v)) dH(u ′ ,v ′ )=d-1 P (u ′ )P (u),
from which (11) is easily deduced. Without loss of generality, one may assume that the states in S are numbered 1, . . . , M and the additional intermediate states in S b \S are numbered M +1, . . . , M b . Then A x d,n (y, z) may be viewed as a matrix indexed by (y, z) ∈ S 2 b and the recursion (11) written more compactly as

Āx d,n+1 = P 0 Āx d,n + P 1 Āx d-1,n , (14)
where Āx d,n = A x d,n (:) is the column vector obtained by stacking the columns of A x d,n and

P h = (P h (z, z ′ )) 1≤z≤M b ,1≤z ′ ≤M b (h = 0, 1), ( 15
) are two M 2 b × M 2 b matrices. Each consists of M b × M b sub- matrices P h (z, z ′ ) = (P h (z, z ′ , y, y ′ )) 1≤y≤M b ,1≤y ′ ≤M b (16)
with

P h (z, z ′ , y, y ′ ) =        0, if z ′ = y ′ (u,v)∈T y ′ y ×T z ′ z dH(u,v)=h P (u), if z ′ = y ′ . (17)
Finally, let I M b be the column vector obtained by stacking the columns of the M b × M b identity matrix. Then (10) may be written as

A d = x∈S p(x) ∞ n=d I T M b Āx d,n = x∈S p(x)I T M b ∞ n=d Āx d,n = x∈S p(x)I T M b Āx d . (18)
V. DIRECT SPECTRUM EVALUATION If the spectral radius ρ(P 0 ) < 1, it is possible to compute the coefficients A d exactly by matrix inversion. This assumption on P 0 is not very restrictive, since if ρ(P 0 ) > 1 there could be an infinite number of non-converged pairs of paths at distance d H = 0, which implies infinite decoding delay. Several JSC-AC FSEs from [START_REF] Ben-Jamaa | Analytical tools for optimizing the error correction performance of arithmetic codes[END_REF] were tested and all had ρ(P 0 ) < 1.

Proposition 3:

The vectors Āx d = ∞ n=d Āx d,n
, leading to the spectral coefficients A d via (18), may be computed as

Āx d = [SP 1 ] d-1 S Āx 1,1 + [SP 1 ] d S Āx 0,1 (d ≥ 1), ( 19
)
where Āx 0,1 and Āx 1,1 are the initial values (12) in vector form, and

S = ∞ k=0 P k 0 , (20)
which converges for ρ(P 0 ) < 1 and may be evaluated as

S = (I -P 0 ) -1 . ( 21
)
Proof: The proof relies on the following. Lemma 4:

Āx d+1 = SP 1 Āx d (d ≥ 1). ( 22
)
To prove ( 22), expand the recursion (14) over k steps in the variable n, from n = d + 1 to n = d + k + 1, and over one step in d, yielding

Āx d+1,d+k+1 = k i=0 P k-i 0 P 1 Āx d,d+i . (23)
In ( 23), the facts that Āx d,n = 0 for n < d and P 0 0 = I have been used. Then expand Āx d+1 as

Āx d+1 = ∞ k=0 Āx d+1,d+1+k = ∞ k=0 k i=0 P k-i 0 P 1 Āx d,d+i = ∞ k=0 P k 0 P 1 Āx d,d + ∞ k=1 P k-1 0 P 1 Āx d,d+1 + ∞ k=2 P k-2 0 P 1 Āx d,d+2 + . . . (24)
The sum (24) may be rearranged into

Āx d+1 = ∞ k=0 P k 0 P 1 ∞ i=0 Āx d,d+i = SP 1 Āx d ,
proving Lemma 4, thanks to which it is sufficient to show

Āx 1 = S Āx 1,1 + [SP 1 ]S Āx 0,1 (25)
in order to prove (19). Using (20), the first term in (25) is recognized as the recursive summation of the first term in (14).

For the second term in (25), first Āx 0 = S Āx 0,1 is obtained in the same fashion as the first term, then the same steps as for Lemma 4 are applied to get [SP 1 ] Āx 0 .

VI. BOUNDING THE APPROXIMATION ERROR

In practice, the spectral coefficients A d are only computed up to a maximal distance d max and, if recursion (11) is used instead of matrix inversion (21), the number of terms in the sum (10) is limited to some N > d max . Hence A d is only approximated from below (because all terms are non-negative) and so also [START_REF] Ben-Jamaa | Analytical tools for optimizing the error correction performance of arithmetic codes[END_REF] can only be approximated. Nevertheless, it is possible to upper-bound the approximation error and thus obtain a proper upper bound. For d ≤ d max and N ≥ d, one may rewrite (18) as

A d = α d (N ) + ε d (N ), where α d (N ) = x p(x)I T M b ᾱx d (N ) = x p(x)I T M b N n=d Āx d,n , (26)
ε d (N ) = x p(x)I T M b εx d (N ) = x p(x)I T M b ∞ n=N +1 Āx d,n (27)
are the N -term approximation and error, respectively. The latter can be bounded with

ε d (N ) ≤ x p(x) εx d (N ) 1 , (28)
where the norm Ā 1 = i | Āi | for vectors Ā and P 1 = max j i |P ij | for matrices P . Using reasoning similar to that leading to (25), one can show

εx d (N ) = SP 0 Āx d,N + SP 1 Āx d-1,N + εx d-1 (N ) , (29)
where the factor P 0 in the first term accounts for the fact that the sum (27) starts with Āx d,N +1 . Because d free dominates the performance at high SNR, a key issue in the recursive approximation of the spectrum via (11) or ( 14) is to determine when d free has been found, i.e., when no unmerged pair of N -bit paths could lower the smallest d such that A d > 0. A sufficient condition was given in [START_REF] Ben-Jamaa | Analytical tools for optimizing the error correction performance of arithmetic codes[END_REF], which adapted to the present framework states that when N is such that A d = 0 and Āx d,N 1 = 0 for d = 0, . . . , d * -1, x ∈ S, while A d * > 0, then d free = d * has been found. Let N satisfy these conditions in the following, implying ε d (N ) = 0 for d < d free . Then (29) may be rewritten as

εx d (N ) = SP 0 Āx d,N + d-1 i=dfree [SP 1 ] d-i S Āx i,N (d free ≤ d ≤ d max ). (30)
For d > d max , the error may be computed using (22) as

εx d (N ) = Āx d = [SP 1 ] d-dmax Āx dmax = [SP 1 ] d-dmax ᾱx dmax (N ) + εx dmax (N ) (31)
Since the norm • 1 is sub-additive and sub-multiplicative, (30) in conjunction with (28) suffices to bound the approximation error for d = d free , . . . , d max using the quantities P 0 1 , P 1 1 , S 1 , Āx d,N 1 and ᾱx dmax (N ) 1 . These are all readily computed during the recursive spectrum evaluation, except for S, which may be upper bounded as follows.

Lemma 5: If the spectral radius of P 0 satisfies ρ(P 0 ) < 1, then for any induced matrix norm • there exists n 0 > 0 such that for all n ≥ n 0 , one has

S ≤ n i=0 P i 0 1 -P n+1 0 . ( 32
)
Proof: Since ρ(P 0 ) < 1, for any induced matrix norm there exists n 0 > 0 such that P n 0 < 1 for any n ≥ n 0 . One may write (20) as

S = n i=0 P i 0 + ∞ i=n+1 P i 0 = n i=0 P i 0 + P n+1 0 ∞ i=0 P i 0 .
Then, for any n ≥ n 0 , the lemma follows from The tail sum for d > d max in the union bound (5) may also be bounded, since the error term (31) grows polynomially, while the channel error probability [START_REF] Pollara | Finite-state codes[END_REF]

γ = E b N0 -log κ. Then ∞ d=dmax+1 A d P d ≤ β ∞ d=dmax+1 κ d-dmax δe -d E b N 0 = βδe -γdmax κ dmax (e γ -1) (γ > 0), (33)
where δ = 1 2 erfc d free

E b N0 exp d free E b N0
stems from tightening the bound (8) as in [START_REF] Viterbi | Principles of Digital Communication and Coding[END_REF]Ch. 4.5]. The tail bound converges only if γ > 0, i.e., if the channel SNR is larger than log κ. In many settings this turns out to be problematic, since SP 1 1 and thus κ may be quite large (κ ≈ 30 for the example FSC in Sec. VII), making it impossible to obtain a useful (tight) bound at practical SNR levels. One may choose to neglect the tail sum if d max is large enough; a slightly preciser method consists in reducing κ to a value that is closer to the asymptotic growth rate κ ∞ = lim d→∞ A d+1 A d . Empirically, this was found to be close to the spectral radius ρ(SP 1 ) already for small values of d > d free + 1 (κ ∞ ≈ 2.6 for the FSC in Section VII).

VII. EXPERIMENTAL RESULTS

Figure 1 shows the encoder graph of a JSC-AC for a Bernoulli source with Pr(X= 0) = 0.11, obtained as outlined in [START_REF] Ben-Jamaa | Analytical tools for optimizing the error correction performance of arithmetic codes[END_REF]. Simulation results (for blocks of 1024 symbols, BPSK over AWGN and ML sequence (Viterbi) decoding) are compared with approximate and proper upper bounds in Figures 2 and3, confirming the empirical "rule of thumb" that relatively few terms are necessary to obtain usable approximate bounds (notice that for large enough N the bound using recursive approximation can be close to that using matrix inversion). Unfortunately, the tail bound (33) provides rather gross overestimates even when using the empirical growth rate κ ∞ .

Fig. 1 .

 1 Fig. 1. Graph of a finite-state joint-source channel arithmetic encoder for a Bernoulli source with Pr(X = 0) = 0.11. Average rate 1.01 bits/symbol, d free = 2, A 2 = 0.44.

Fig. 2 .Fig. 3 .

 23 Fig.2. Recursive approximation of the union bound: partial sum of (5) up to dmax = 10, using partial sums of (10) up to n = N , where N = 10, 20, . . . , 100 (evaluated with (11)); proper upper bound including bounded error terms (30) and tail sum (33) with dmax = 10, N = 100, κ∞ = 2.6.

There is little practical interest in FSEs with multiple ergodic components, however, all results apply if appropriately averaged over these components. Transient components may be eliminated, since they do not affect steady-state performance.

The tools developed here can also be applied to periodic codes, either by carrying out all computations for an entire period (bn, . . . , b n+m-1 ) at a time, instead of a single bit bn, or by separate computations for each of the m phases.

NEWCOM++ (contract n. 216715


Comment on "Large Slip of Aqueous Liquid Flow over a Nanoengineered Superhydrophobic Surface" In a recent Letter [1], Choi and Kim reported slip lengths of a few tens of microns for water on nanoengineered superhydrophobic surfaces, on the basis of rheometry (cone-and-plate) measurements. We show that the experimental uncertainty in the experiment of Ref. [1], expressed in term of slip lengths, lies in the range 20 -100 micrometers, which is precisely the order of magnitude of the reported slip lengths. Moreover we point out a systematic bias expected on the superhydrophobic surfaces. We thus infer that it is not possible to draw out any conclusion concerning the existence of huge slip lengths in the system studied by Choi and Kim.

Choi and Kim performed torque measurements using a commercial rheometer (AR2000, TA Instruments) with a cone-and-plate geometry. The slip length δ is deduced from the correction of the torque M w.r.t. a reference noslip value

M 0 : δ = 2 3 θ 0 R(1 -M/M 0 ). M 0 = 2π 3 µΩR 3 /
θ 0 is the prediction in the absence of slippage; R is the cone radius (≃ 3 cm) and θ 0 its opening angle (2 • ); µ is the viscosity of the liquid; Ω is the rotational velocity. Putting numbers in this expression shows that the reported 20µm slip length for water corresponds to a 3% correction to the (small) reference torque M 0 ∼ 5µNm. The authors however claim a 3µm uncertainty on the slip length, which corresponds to an overall 0.5% errorbar on the relative deviation of the torque (M 0 -M )/M 0 . Such a precision is not attainable in the present experiment. To illustrate this uncertainty issue, we have performed benchmark experiments using an AR2000 rheometer with a smooth, stainless steel cone-and-plate geometry (with the same radius and cone angle), very close to that of Ref. [1]. This rheometer was calibrated using a reference Newtonian silicon oil (BR0050CPS, µ = 48.4 mPa s at 25 • C), which yielded the expected value to within 0.5%. Then, turning to distilled deionized water at 25 • C, we performed torque measurements similar to [1] and measured the torque standard deviation ∆M/M with this liquid, yielding ∆M/M ≃ 1.4% at γ = 150 s -1 , up to 3.5% at γ = 50 s -1 . This is far above the claimed 0.5% uncertainty. If the uncertainties on the filling volume and on the gap size are included, the global uncertainty on the measured torque is at least ∆M/M ≃ 3%. Finally, using the expression δ(M ) and adding a similar uncertainty on the reference viscosity µ in M 0 , leads to ∆δ/δ ≃ 100-200%, so that ∆δ ∼ 20-40µm for water and ∆δ ∼ 50 -100µm for glycerin. The reported effect are therefore within uncertainty and the experiment of Ref. [1] should be considered as inconclusive.

Another source of difficulties in the interpretation of the experiment in Ref [1] is the role of secondary flows. The relevant reynolds number is Re = ρΩR 2 θ 2 0 /µ 2, and inertial effects should lead to an increase in the torque up to two percent (M -M 0 = 6.10 -4 Re 2 [START_REF] Walters | Rheometry[END_REF]). In view of the claimed resolution, this correction should be measurable in the experiment of Ref. [1] and interpreted as an apparent negative slip length up to ∼ -15µm for the smooth hydrophilic surfaces. This effect is not detected in [1], which furthermore confirms the weakness in the interpretation of the measurements. At this stage, it is worth pointing out a systematic bias on the superhydrophobic surfaces. Indeed, for the same liquid volume filling the gap in the cone-and-plate, the meniscus at the edges makes the radius R slightly smaller on the super-hydrophobic surface (with very large contact angle), w.r.t. the other surfaces with smaller contact angles. The variation in R is predicted to be of the order of a fraction of the gap at the border ∆R ∼ -αθ 0 R (with α 1). Assuming no-slip at the surfaces, the resulting decrease of the solid-liquid area leads to a reduction of torque on the superhydrophobic surfaces, ∆M cap /M 0 ∼ -3αθ 0 , of the order of a few percents. The misinterpretation of this effect using Eq. ( 1) thus erroneously predicts a slippage effect with a slip length δ eff ∼ αθ 2 0 R, of the order of a few tens of micrometers.

In summary, the experimental uncertainty that we estimate is comparable to the amplitude of the effect the authors have observed. Moreover a systematic bias could be wrongly interpreted in terms of very large slippage on superhydrophobic surfaces. The experiments of Ref. [1] are therefore inconclusive.

  Lydéric Bocquet 1 , Patrick Tabeling 2 , Sébastien Manneville 3

	(1) LPMCN, Université Lyon I,
	UMR CNRS 5586, Lyon, France
	(2) PCT, ESPCI, UMR CNRS 7083, Paris, France
	(3) CRPP, UPR CNRS 8641, Pessac, France


Pierre Seppecher

MOVING CONTACT LINES IN THE CAHN-HILLIARD THEORY

We establish the equations of motion of an isothermal viscous Cahn-Hilliard fluid and we investigate the dynamics of fluids having moving contact lines under this theory. The force singularity arising in the classical model of capillarity is no longer present. This removal is due to a mass transfer across the interface combined with a finite thickness of the interface. A numerical simulation of the flow in the immediate vicinity of the contact line shows the connection between the static contact angle, the dynamic angle and points out the influence of the velocity.

INTRODUCTION

Capillary phenomena are important in many circumstances, for instance, the shape of the interface of a drop lying on a smooth horizontal surface depends on gravity and surface tension, but to compute 1:his shape, one needs to know another parameter: a boundary condition for the interface must be written. In general, the angle made by the interface and the surface is given. This contact angle is the parameter which models the physical interactions between the surface and the fluids. These interactions may have great influence upon the statics or the dynamics of the fluids. For instance, they are responsible for the rising of a fluid in a capillary tube. When the interface is moving with respect to the surface, these interactions cannot be modelled in a simple way: it is well known [1][2][3] that the assumptions that the phases are both viscous fluids and that the no-slip condition holds on the boundary lead to a force singularity at the contact line. The fundamental reason for this paradox is that the no-slip condition combined with the balance of mass across the interface implies a discontinuity of the velocity field at the contact line: the velocity cannot belong to the functional space H ~. The total dissipation which is equivalent to the H~-norm is infinite. It must be noted that the problem remains when allowing mass transfer across the interface [START_REF] Seppecher | Congrds Frant:ais de M~canique[END_REF]. Another problem must be emphasized: the curvature of the interface is related to the pressure jump and, according to the theory, this jump tends to infinity when approaching the contact line. Such a situation is physically impossible.

There are different ways to overcome this difficulty. A first one is to consider that a thin film made by one of the phases covers the surface so that the interface does not reach the surface [START_REF] De Gennes | [END_REF]. Then the apparent contact line can be studied using classical continuum mechanics but the difficulty arises again at the tip of the film.

The second one, the most popular, is to relax the no-slip condition by introducing a new parameter, the slip length [6][7][8][9]. This parameter as well as the dynamic contact angle must be determined experimentally. They may depend on the velocity of the contact line. The difficulty is to make a distinction between the contact angle and an apparent contact angle [10,11].

The third one, frequently used, is to assume that one of the phases is a perfect fluid. The paradox disappears. The flow in the viscous phase near the contact line is a rolling-like motion, but however small the viscosity coefficients of the phases may be, the dissipation is infinite. Assuming that one of these coefficients vanishes means passing carelessly to the limit in a singular perturbation problem. Some important phenomena may be missed.

Finally, the continuum model used may be questioned. It is itself an approximation which may not be valid in the vicinity of the contact line. Then one may use molecular simulations 978 P. SEPPECHER [12] or a continuum model able to describe the interface as a layer of finite thickness. This is the way chosen in this paper. We use the Cahn-Hilliard model (or the van der Waals model). This model was not built in order to solve the contact line paradox; its parameters can be estimated from situations irrelevant to the problem. It is clear that it is a simple model probably superseded in most practical situations but it is the first continuum model which describes the flow at a moving contact line.

The statics of the Cahn-Hilliard fluid is well known, but the nature of the internal forces in such a fluid is not trivial. One needs the second gradient theory [13] or the theory of continua with edge forces [START_REF] Dell Isola | Mecanica[END_REF] to understand it. In the following section we establish the equations and boundary conditions for a viscous isothermal Cahn-Hilliard fluid. In Section 3 we restrict our study to the vicinity of a contact line. We define the zone (the inner zone) where the Cahn-Hilliard model will be used. In an intermediate zone where the classical model can be used we develop an analytic solution. Then in Section 4 we show how this intermediate solution matches the external and inner solutions and fixes the boundary conditions for the inner problem. We point out the dimensionless parameters for the inner problem. This problem is non-linear. We study it by numerical simulation which is the subject of Section 5, where we split the problem into two minimization problems. The minimization of the dissipation is a linear problem solved in a classic way while the minimization of the energy is a non-linear problem solved by a steepest-descent method. The dependence of the dynamic contact angle upon the velocity of the contact line is clearly exhibited.

THE ISOTHERMAL VISCOUS CAHN-HILLIARD FLUID

Free energy, dissipation and constitutive laws

The model of Cahn and Hilliard [START_REF] Cahn | [END_REF] is the simplest continuum model for multi-phase fluids. It is assumed that the volume free energy E of a fluid whose mass density is denoted p is the sum of a non-convex volume energy W(p) and a term taking into account the non-homogeneity of the fluid

W(p) + ~ (Vp) 2. (2.

1) E

The function W is a two-well positive function vanishing only for p e {pv, Pt}. Because of the shape of W, the fluid tends to divide into two phases p = Pv and p = pt and the term (,~/2) ($,p)2 tends to reduce the variation of the field p, turning the interface into a thin layer and endowing it with an energy (the surface tension).

The most interesting theoretical feature of this model is the fact that it is not compatible with classical continuum mechanics [16]. A direct way to point out this fact is to write the Clausius-Duhem inequality without making any assumption upon internal forces power p~n,, energy flux J~ or entropy flux J~. Under isothermal conditions (T = constant) this inequality reads where : and i denote the contraction products of tensors of order two and three.

It is then clear that inequality (2.2) is incompatible with the assumptions that the entropy flux is equal to the energy flux divided by the temperature and that the power of internal forces is the product of Cauchy stress tensor times velocity gradient. This incompatibility can be removed by changing: the form of the energy flux [16]; or the form of the entropy flux following MUller [START_REF] Misller | Thermodynamics[END_REF]; or t]~e form of the power of internal forces [START_REF] Casal | [END_REF][19][20].

These differerLt ways are formally equivalent but we claim that the last is the most natural one as, for a Cahn-Hilliard fluid, internal forces cannot be represented by the Cauchy stress tensor alone. Indeed a fundamental postulate in Cauchy's construction of the stress tensor (using a tetraheclron) is that no line forces are present on the edges of a domain: this postulate is violated by a Cahn-Hilliard fluid [21]. On the other hand, when studying the equilibrium of a Cahn-Hilliard fluid, i.e. the minimization of free energy, an extra boundary condition appears which cannot be explained by means of the Cauchy stress tensor.

Then we assume that J~ = TJs and that pint has the form pint = _ T : grV -C i VVV.

The Clausius-Duhem inequality becomes F = (r-r"):VV + (C -C"):VVV-> 0.

(2.5)

In the approximation of linear thermodynamics the dissipation F is a non-negative quadratic form of (VV, V~V). Even with this approximation the general constitutive laws involve many coefficients so we admit that there is no dissipation due to VVV and no anisotropy due to the presence of the thermostatic parameter Vp. Then the constitutive law involves only two viscosity coefficients, v and ~, and the dissipation has the usual form

F = v tr(D) 2 + 2~D: D, (2.6)
where D = ½(VV + VV t) and stress tensors r and C verify the constitutive laws: C = C", r = r" + 11 with II = v tr(D)Dd + 2~D.

(2.7)

Balance of forces and boundary conditions

The easiest way to write the balance of forces and the boundary conditions for a Cahn-Hilliard fluid is to use the virtual power principle in the framework of the second gradient theory [13]. An alternative (and equivalent) way, closer to Cauchy's approach, is to consider the fluid as a continuum in which edge forces are present [START_REF] Dell Isola | Mecanica[END_REF]. This makes the nature of the boundary conditions clearer. The external medium (the container) exerts on the fluid three types of contact force distributions (in this paper we do not need to consider any body force): a surface density of forces on the whole boundary; a line density of forces on the edges of the boundary (if any); and a surface density of double forces (i.e. a distribution of order one with respect to transverse derivatives) on the boundary (see [START_REF] Dell Isola | Mecanica[END_REF] for a discussion on the nature of this distribution).

We admit that the no-slip condition holds on the whole boundary. Then, when the boundary is not moving, the velocity of the fluid V vanishes on the boundary. The power of contact force distribution ~xt is reduced to the power of double forces and takes the form (

2.14)

The last equation shows that a Cahn-Hilliard fluid can only support double forces which are perpendicular to the boundary. This is due to the fact the free energy depends only on p (it is a model of fluid) and to our assumption that the dissipation does not depend on VVV. This phenomenon is analogous to the fact that a perfect fluid cannot support contact forces which are not perpendicular to the boundary. In a sense our model is a model of a "semi-perfect" (2.15) Before using the system of equations (2.13), (2.15), let us see its connection with the classical theory of capillary phenomena.

Connection with the classical theory of capillarity

The connection between the model of Cahn-Hilliard and the classical model of capillarity has been established rigorously in the static case [22].

The model has its own characteristic length L [see definition (4.9) of L]. This length is actually characteristic of the thickness of the interface. It is in general much smaller than the characteristic size of the container: a small dimensionless parameter is introduced. Then an asymptotic study of the model of Cahn-Hilliard is possible. This study was carried out first by Cahn and Hilliard [23], and rigorously, using the notion of F-convergence, by Modica [22]. This study has only been made in the static case where the problem is to find a minimizer p of the free energy. The boundary condition (2.15) is taken into account by an extra energy term:

-f;,c~ AGp ds. Some standard assumptions on the behaviour of all physical quantities with respect to the small parameter (see [START_REF] Bouchitte | Motion by Mean Curvature (Proceedings)[END_REF][START_REF] Seppecher | [END_REF][26] for examples where these assumptions do not hold) are necessary to find the expected limit model which can be summarized by: (2.18)

tr:~ = --A Gp4 + dp if G > 0, (2.19) = -AGp3 + I '°' 2X/~--W(p) dp if G < 0. (2.20) or.~ "tO3
The static contact angle ~s is then determined by the Young law: When equation (2.16) does not have four solutions or when equation (2.21) has no solution, the wetting is complete. We will not consider these cases in this paper.

To our knowledge no result has been obtained in the dynamical case. We think that this asymptotic model is likely to be valid far from a moving contact line. This paper is devoted to the study of the vicinity of the contact line.

DESCRIPTION OF A MOVING CONTACT LINE

line. On these two assumptions we can study the flow near the moving contact line without considering the flow in the whole domain.

While the flow in the external zone is three-dimensional, we assume that it is twodimensional in the inner and intermediate zones (here we assume implicitly that the curvature radius of the contact line is much larger than the characteristic size of the inner or intermediate zones), so we reduce our study to a plane perpendicular to the contact line (the plane of Fig. 1). Moreover, we assume that, in a system of coordinates tied to the contact line, the interface is not moving and the flow is stationary. In such a system of coordinates it is the boundary of the domain which is moving (with a velocity tangent to the boundary). Only the normal component (normal to the contact line) of this velocity is well defined: we can state that this velocity lies in the plane perpendicular to the contact line. We denote it Vcx~.

More precisely, we assume that (in a vicinity of the contact line) every tangent line (O, P) to the interface intersects the boundary. Let R~ be a length much larger than the characteristic thickness of the interface and much smaller than the characteristic length of the external flow, and let Re be a length much larger than R~ but much smaller than the characteristic length of the external flow. Let us consider the tangent line (O, P) such that IIOPII = R~ and let us choose O as the origin of the system of coordinates (O, Xl, x2) (x~ is a unit vector tangent to the boundary and normal to the contact line, see Fig. 1). This system of coordinates is tied to the contact line. Using polar coordinates (r, 0), we define the inner zone as the set {r<R~, 0 < 0 < 7r} and the intermediate zone as the set {R~ < r < R2, 0 < 0 < if}. The angle (xt, OP), denoted qb, is called the apparent contact angle.

The curvature radius of the interface in the external zone is much larger than R2. We assume that this is also true in the intermediate zone, thus the interface is assumed to be near plane in this zone. In fact R~ must be large enough because we know that the curvature increases (diverges) in the classical theory of capillarity when approaching the contact line.

The following section is devoted to the study of the flow in the intermediate zone in order to obtain a boundary condition for the inner zone.

Intermediate flow

In the intermediate zone we use the classical theory of capillarity. The interface is a near plane surface and no mass transfer occurs. The two phases are homogeneous and incompressible. For the sake of simplicity we assume that their viscosity coefficients are identical (no particular difficulty should arise if different viscosity coefficients are considered).

In the plane (O, xl,x2) we use polar coordinates (r, 0) and we study a stationary two-dimensional flow in the domain @ = {0 < 0 < 7r, R i < r < R2}. We assume that two phases M = {0 > ~} and ~ = {0 < dO} with mass densities p.~ and p~ are divided by a plane interface # = {0 = ~} (see Fig. 2). The normal to # pointing into M is denoted n.

Let V.~.~, V~, r.~, and r~ be the velocities and stress tensors in each phase. Let us recall our assumptions:

(i) both phases are homogeneous incompressible fluids:

(r~" n), = (r:~. n), on #, (

where II denotes the projection to the plane tangent to the interface; This solution cannot be the physical solution for a moving interface in the vicinity of the contact line (R~ = 0) as the dissipation in any open part of the domain containing O would be infinite, and the normal stress jump across the interface would diverge when approcaching O, thus the Laplace', law binding this jump to the interface curvature could not be satisfied even approximately. However, this solution can be a good approximation for the flow in an intermediate zone (far enough from O for the normal stress jump to be small and close enough to O for the external flow to have a weak influence). We will use it in Section 4 to write the boundary conditions for the inner problem. For this purpose, we denote by Vv~'°~"°~'*(0) the velocity on ~,~ t.J ~ corresponding to W~ or W~ depending on whether 0 < • or 0 > • (the superscripts recall the dependence of this boundary velocity field on the parameters V~, D~, D:~ and 0).

INNER PROBLEM

In the inner zone the characteristic length is so small that the classical theory of capillarity is no longer valid: we cannot neglect the thickness of the interface. A non-zero thickness is fundamental because it removes the discontinuity of the mass density field, then it removes the discontinuity of the velocity field when mass transfer occurs, and finally removes the singularity of the dissipation at the contact line.

Boundary conditions

The inner domain and the expected solution are represented in Fig. 3. The domain is a half disc @n, limited by a half circle Yn, and a segment 5¢R,. boundary. In this drawing the thickness of the interface (i.e. the zone where the isodensity curves are concentrated) was exaggerated. In fact the interface at r = R~ looks like a line in an R~-scale. Let us discuss the static boundary conditions first: (i) Due to our definition of inner zone (see Fig. 1), the interface intersects ZR, orthogonally. It is known from the equilibrium theory of a Cahn-Hilliard fluid ]equation (2.22)] that this condition is obtained by setting G = 0 on ZR,.

(ii) At equilibrium the contact angle ~ is connected with the value of G on 5eR, ]equation (2.21)]. This value is given by the physics of wall-fluid interactions. Fixing G is equivalent in the Cahn-Hilliard theory to fixing the static contact angle in the classical theory of capillarity.

(iii) At this point the static problem is not well-posed. We must write an extra condition deriving from the external solution. We have already discussed this point in Section 3.1 and concluded that the interface curvature should be fixed. As such a condition raises technical difficulties we prefer to fix the total amount M of mass lying inside the inner zone. The parameter M will then be adjusted to make the curvature of the interface vanish at r = R~ (therefore the interface fits with a plane interface in the intermediate zone).

Conditions (i)-(iii) lead to a well-posed static problem whose solution corresponds to a plane interface intersecting the wall at the origin with an angle compatible with Young's law. This situation is used for a first validation of the numerical simulation we describe in Section 5. It must be noted that, due to its non-linearity, this problem in general does not have a unique solution.

Let us now consider dynamical conditions: (i) On b°R, we admit the no-slip condition V = Vcx~.

(ii) On the boundary ZR,, V coincides with a solution V v''D.'°*'° obtained for the intermediate zone (see Section 3.2).

As the parameters ~, D, and D:~ are unknown (the study of • is one of our main objectives), we must write some supplementary conditions.

(iii) We deal only with mass density fields p whose variations are concentrated in a thin layer. Then we can measure the apparent contact angle (see Fig. 3 or Fig. 4) and we state that qb coincides with this angle.

(iv) As we study stationary flows, the total mass supply in @R, must vanish (fx~, P Vv~'°'D~'*" n dl = 0), so:

D,=-[fx~, pVV~'°'t'*'n d/I-' [fx~, pvV~'L°'*" ndl]D,~. (4.1) (v)
The mass transfer from phase M to phase ~ is unknown. This uncertainty comes from external conditions (see Section 3.1). We assume that the phases are close enough to the chemical equilibrium for the mass transfer across the interface to be negligible for r-> R~. Thus D.~ will be adjusted to make mass transfer vanish at r = Rt. These last three conditions allow us to compute the parameters ~, D~ and D.,~ when the field p is known. Then, using (ii), we know the velocity field on the boundary Zn,.

Let us sum up our boundary conditions: we write a non-homogeneous Neumann condition for the mass density field p, plus a global condition upon the mass M contained in the inner zone. For the velocity field, we write a non-homogeneous Dirichlet condition which depends on the mass density field and on a mass transfer D.~¢. The two parameters M and D~ represent the influence of external conditions upon the flow in the inner zone. An adequate selection of these parameters will lead to a near plane interface without mass transfer at r = R~.

Equations of the inner problem, dimensionless parameters

Recalling the set of equations stated in Section 2 and the boundary conditions discussed in the previous section, the problem is to find a vector field V and a scalar field p such that: pV( oW-AAp] = V. (II) in DR, (vi) finally we adjust parameters m and d u and repeat the whole procedure in order to have a near plane interface without mass transfer when r = R,.

(U q-Um)V( ~--Au) =GAY" (n) V'((U+Um)V)=0 inDR
Problems (5.1) and (5.6) are discretized by a finite difference method. Problem (5.10) is linear. We solve it by using a bi-conjugate gradient method. The problem (5.1) is strongly non-linear. We solve it by a steepest-descent method. We know that, in general, it does not have a unique solution. In order to find the expected solution we need to initialize the mass density field u with a field which is not too far from the solution. This is not very difficult for the equilibrium configuration and, in the dynamical case, we use as an initializing field the solution u obtained for a slightly smaller capillary number Ca. Otherwise the most likely disaster is the inversion of the position of phases M and ~ which tends to happen when one of the phases tends to wet the boundary completely.

Dependence of the dynamic contact angle on the parameters

The dimensionless free energy w is a positive function such that w(-1)= w(1)=0 and fl_~ ~(t) dt = 1. In our simulation we use the polynomial w(u) = ~(1 -u2) 2. Such a form for the free energy is valid when the fluid is not too far from the critical point. We hope that, qualitatively, our results do not depend on the form of w.

In Fig. 4 we show the density field in a typical situation: R = 20, Um----2, K = 10, g = 0, Ca--20 × 10 -3. In Fig. 5 we show the flow lines in the same situation. We notice that the curvature of the interface and the transfer of mass across the interface are concentrated in the immediate vicinity of the contact line, so the dynamic contact angle is well determined. Here it is about 1.9 while the static contact angle is Ir/2.

In Fig. 6 we show the dependence of the contact angle on the capillary number Ca in two cases (R=20, um=2, K=10, g=0andR=20, Um=2, K=10, g=-0.3). When Ca=0the apparent contact angle coincides with the theoretical static angle predicted by equation (2.21), then its increase is almost linear with respect to the capillary number. The growth rate does not depend on the static angle (see Fig. 6).

We have noted a weak influence of parameter urn: the growth rate decreases slightly when Um increases. This dependence can be interpreted in the following way: when approaching the critical temperature, Um increases, the difference between the mass densities of the phases decreases and the mass transfer is easier.

As expected, the growth rate depends strongly on the parameter K which is responsible for ,Fig. 4. Density field in the inner zone. mass transfer dissipation. Figure 7 shows the influence of K in the situation of R = 20, /4 m ~-2, g=0, Ca=20×10 -3.

CONCLUSION

The Cahn-Hilliard model used in this paper is very simple and comparisons with experiments would be difficult. However, it gives an enlightening view of the flow at a contact line. First of all, mass transfer is the fundamental fact which removes the force singularity (in this connection it would be interesting to extend this study to the motion of an interface dividing two slightly mixable fluids). Secondly, it shows that the main change in the direction of the interface lies in the immediate vicinity of the contact line (the same region where mass transfer occurs). In this region the interface cannot in any way be considered a dividing surface and the classical model of capillarity cannot be applied. Thus, any study considering the interface a bidimensional surface must use a velocity-dependent contact angle: our results are in favour of model "D" following the notation and discussion of Hocking [11]. Indeed, no real angle equal to the static angle can be defined. We found a near linear relation between the angle and the velocity of the contact line. This linear dependence gives reasonable results if there is no hysteresis of the static contact angle [8]. Indeed, there is no hysteresis in our model. The way to model contact angle hysteresis in the Cahn-Hilliard theory is to consider rough surfaces. This has not been done here (and would be incompatible with the steady flow assumption).

For a large capillary number, some problems of convergence of our algorithm arise. That is why in Fig. 6 the apparent contact angle does not reach the values ~ = 0 or • = Jr. It might be a purely numerical phenomenon. It might be connected with the non-uniqueness of the solution of (5.1); but it might also be due to the physics: a steady flow may not exist for a large velocity or a small angle. Indeed experiments (as reported in [9]) seem to point out a maximum value for the contact ~Lngle.

Our study is valid in isothermal conditions. It is also valid in adiabatic conditions (it is enough to replace the free energy E by the internal energy); but the energy transport may be the factor which limits phase transition near the moving contact line. It should be necessary to add the heat equation to our problem (4.3)-(4.8).

(2. 2 )

 2 F=T~7. J~-J~ -Pint-p dtLet us define the following tensors of order two and three: r ''= W(p)-p~pthe identity tensor. Using the balance of (dp/dt = -pV. (V)) we can compute d(E/p)/dt and we get: d(E/p) p ---~" : VV + C" i VVV (2.4) dt

H 12 )

 12 ~xt= where n denotes the external normal to the boundary O@. Let ~/be the acceleration of the fluid; the principle of virtual powers reads f pv.vdv=~:.:H.(n. Vv)ds-f (r:Vv+CiVVv)dv (2.9) for all fields v such that v = 0 on ~99. Using the divergence theorem, this equation becomes f[ov -v. (r -v. (c))] • v dv -f~, [n -n. C. n]. (n. Vv) ds = Taking the constitutive laws (2.7) into account, the equations of motion of the isothermal viscous Cahn-Hilliard fluid read p~l = -pV(O~p-AAp)+ V.(H) in ~, (2.13) H= -Ap(n. Vp)n on 09.

  fluid. Defining G by H = -ApGn, equation (2.14) reads n. Vp = G on O@.

  (i) the fluid is divided into two incompressible homogeneous phases ~ and ~, whose mass densities are pv and pt; (ii) a surface energy tr = f~'v ~ dp is associated with the surface Lr dividing ~ and (the interface); (iii) surface energies tr~e and tr:~ are associated with the boundaries 09 71 ~/and c99 71 ~ of the phases. Let pl </32 <~ p3 ~ p4 be the four solutions of the equation X/2A-tW(p) = IGI. (2.16) Then o'.u and tr:~ are given by o'.e = -AGoz + 2V~-W(0) d 0 if G >0, (2.17) v tr~¢ = -AGpl + dp if G <0.

  cos(~) = tr.~ -try (2.21) or Example When G = 0 'we have p~ = p2 = pv and P3 = P4 ----Pl, SO O'~ = O'.o~ ----0, hence ~s = ~. (2.22)

Fig. 1 .

 1 Fig. 1. Partition of the domain.

Fig. 2 .

 2 Fig. 2. Flow in the intermediate zone.

  xdo sin(do) 2 a.u = V~do ---~-)~)(-~-2~-)--si~)' (3.16) -(to -do) -sin(do)cos(do) a:~ = Vc do(l r _ do) _ sin(do)cos(qb)(tc -2do) -sin2(do) ' -ffdo sin(do)cos(do) + sin(do) 2 b.~ = V~ do(x -do) -sin(do)cos(do)(zr -2do) -sin2(do) ' (3.17) (3.18) -do sin(do) 2 c.~¢ = V~ do(tr -do) -sin(do)cos(do)(x -2do) -sin2(do) ' (3.19) (tr -do)sin(qb) 2 c:, --V¢ do(zr -do) -sin(do)cos(do)0r -2do) -sinZ(do) ' -dp cos(do)sin(do) + sin(do) 2 d,¢ = V~ do(zr -do) -sin(do)cos(do)(ff -2do) -sinZ(do) ' (x -do)cos(do)sin(do) + sin(do) 2 d:~ V~ do(jr -dO) -sin(q~)cos(do)0r -2do) -sinZ(do) ' 3(27r -do) -3 g.~ = u.~2---7--, _-L~rUr ~-z 4-D,~--(3.23) ' 27rdo' 3 -3(,I, + x) g:~ = D~ 27rQr -do) + D~ 2zrdo2 , (3.24) -3(~r -O) -tr 3 h~ = D,,, 21r(Tr -0) 3 + D,~ 2trO(~r -O)= D.~2xO(zr _ O) + D.~ 2zrO---------5-. (3.26)

Fig. 3 .

 3 Fig. 3. Geometry of the inner problem, expected solution.

11 )

 11 , M, D~ and D~ being such that (i) • is the apparent contact angle; (ii) the curvature vanishes at r = R~; (iii) no mass transfer occurs at r = R1; (iv) total mass supply vanishes.We define the characteristic surface energy or, mass density Pd and length L is characteristic of the thickness of the interface and o-is the surface tension. We define the By construction R >> 1. Far enough from the critical point u m is of order 1. g is of order 1 and fixes the static contact angle. The capillary number Ca is the ratio of the viscous forces to the surface tension. It is the most important parameter.Using the dimensionless mass density u = (2pd)-l(2p-(Pl + Pv)), the dimensionless velocity v = V~-1V and 1Lhe dimensionless free energy w(u)= pZAo-ZW(p), the problem (4.3)-(4.8) reads

  OnZR v = v]'d~"d~'*(0) on Zn f uds=m R where II = K tr(D)Dd + 2D and D = ½(Vv + 7v').

Fig. 5 .

 5 Fig. 5. Flow lines in the inner zone.

Fig. 6 .Fig. 7 .

 67 Fig.6. Dependence of the contact angle on the capillary number.

  991

3.1 Geometricai' description, localizationWe consider a contact line moving on a plane rigid surface and we divide the domain into three parts. The first part, called the "external zone", is the region far away from the contact line where the whole geometry of the domain and the external forces determine the flow. The classical model of capillarity can be used there. We are not interested in the flow in this region we assume that it has a weak influence upon the flow near the contact line. The second part, called the "inner zone", is the vicinity of the contact line; its size is so small that the thickness of the interface cannot be neglected there, the Cahn-Hilliard model has to be applied. We assume that there exists a third region called the "intermediate zone" close enough to the contact line for the influence of the external flow to be weak and far enough from the contact line for the classical model of capillarity to be valid. However, it is clear that external conditions may have an influence upon the behaviour of the whole interface and therefore upon the flow in the inner zone. Let us consider, for instance, the equilibrium of a small bubble lying on a plane. For a given static contact angle, the curvature of the interface depends on the size of the bubble: the shape of the interface in the inner zone depends on it. In order to get rid of this influence we assume that the curvature radius of the interface in the external zone is much larger than the characteristic size of the inner or intermediate zones.On the other hand, in the dynamical case, a possible jump of chemical potential between the two phases in the external zone leads to a transfer of mass across the interface, i.e. a non-vanishing normal velocity on the interface. In order to get rid of this influence we assume that this normal velocity in the external zone is much smaller than the velocity of the contact

Acknowledgements--This research has been supported by the Centre National d'Etudes Spatiales (Contract No 92/CNES/0295) and by the Region Provence-Alpes-C6te d'Azur.

NUMERICAL SIMULATION

Two auxiliary minimization problems

Let p be a given scalar field in DR, and consider the following minimization problem:

-g u u (5.1) inf w(u) -pu

n. Vu -g = 0 on 6en, n. Vu = 0 on YR.

(5.4)

By differentiation (5.3) implies:

V(~uu(U)-Au)=Vp in@R.

(5.5)

Now let u be a given scalar field, let qb e ]0, 7r[, de and d~ be three given parameters such that fx, (u~ + u)Vl "a'~'a.`'*. n dl = 0 and consider the following minimization problem: V" (H) -(u + Um)V q = 0 (5.7)

7" ((Urn + U)V) = 0 in ~R (5.8) v = xl on 5eR (5.9) v = Vl ''c'a."'* on ZR.

(5.10)

Here the field q is the Lagrange multiplier of the constraint V •((Um + u)v) = 0. It is defined up to a constant. Setting p = C,q, equation (5.7) may be written as

Car" (II) --(U + Um)V p = 0.

(5.11)

Algorithm

It is clear from (5.8)-(5.11) and (5.4)-(5.5) that solving simultaneously problems (5.1) and (5.6) leads to a solution of the inner problem (4.12)- (4.18). Such is the principle of our algorithm:

For given parameters m and d~: (i) we first set p = 0, solve (5.1) and get an equilibrium configuration u;

(ii) we compute the apparent contact angle qb, the value of d.~ such that the total mass supply vanishes and then we get the boundary velocity field V~"l"'l"'~';

(iii) we solve (5.8)-(5.11) and get a new field p;

(iv) we solve (5.1) and get a new mass density field u;

(v) then we repeat steps (ii)-(iv) until they converge;


C Pideri

P Seppecher

A second gradient material resulting from the homogenization of an heterogeneous linear elastic medium

Homogenization may change fundamentally the constitutive laws of materials. We show how a heterogeneous Cauchy continuum may lead to a non Cauchy continuum. We study the effective properties of a linear elastic medium reinforced periodically with thin parallel fibers made up of a much stronger linear elastic medium and we prove that, when the Lamé coefficients in the fibers and the radius of the fibers have appropriate order of magnitude, the effective material is a second gradient material, i.e. a material whose energy depends on the second gradient of the displacement.

Introduction

Continuum mechanics is usually understood as a homogenized description of materials which are heterogeneous at the microscopic level. Then, it is natural to expect from any general theory of continuum mechanics to be stable by homogenization procedures. We prove in this paper that the class of Cauchy continua does not enjoy this stability property. Indeed, we show that the effective properties of some periodic elastic material have to be described by a second gradient theory.

We consider a composite material made up of an elastic matrix reinforced with elastic fibers. Both materials are isotropic linear elastic materials, the Lamé coefficients in the fibers being larger than in the matrix. The structure is periodic: we assume that the fibers are parallel cylinders with the circular section arranged along a square lattice (see Fig. 1).

Homogenization procedure consists in studying the limit behaviour of the material when the period of the structure tends to zero. What is the behaviour of the other physical quantities as the period tends to zero? The effective properties of the material strongly depend on them: when the elasticity coefficients in the fibers are of the same order of magnitude as in the matrix and when the radius of the fibers is of the same order of magnitude as the period, the problem is a classic one in homogenization theory: the effective material is still a linear elastic material whose coefficients can be expressed in terms of the geometry and of the elasticity coefficients of the matrix and the fibers [START_REF] Sanchez-Palencia | Non Homogeneous Media and Vibration Theory[END_REF]. We study a different case: we want to describe a composite medium reinforced by very thin and very rigid fibers. Then, it is natural to assume that the radius of the fibers tends to zero faster than the period and that the elasticity coefficients in the fibers tend to infinity.

Let us now fix some notations: by convention, we choose the characteristic length of the domain as the unit length. The period of the lattice is denoted by ε. We study the limit ε → 0 and every quantity which is not assumed to be constant as ε tends to zero, is indiced by ε. For instance, the radius of the fibers is denoted by r ε , the Lamé coefficients in the fibers are denoted by λ ε and µ ε while the Lamé coefficients in the matrix are denoted by λ 0 and µ 0 . Then our assumptions read

r ε ε → 0 , λ ε → ∞ , µ ε → ∞ .
This situation has already been studied by D. Caillerie [START_REF] Cailllerie | Homogénéistation d'un corps élastique renforcé par des fibres minces de grande rigidité et réparties périodiquement[END_REF] who, setting λ ε = (r ε /ε) -θ , µ ε = (r ε /ε) -θ , considered in two cases the limit (ε, r ε /ε) → (0, 0): (r ε /ε → 0 then ε → 0) and (ε → 0 then r ε /ε → 0). He found that both cases lead to an elastic material but that the homogenized elasticity coefficients depend on the limit procedure: the two limits ε → 0 and (r ε /ε) → 0 do not commute. Here we let r ε ε , µ ε -1 and λ -1 ε tend to zero together and assume that:

lim ε-→0 r ε ε = 0, lim ε-→0 ε 2 log(r ε ) = 0, lim ε-→0 µ ε r 4 ε ε 2 = µ 1 > 0, lim ε-→0 λ ε µ ε = .
This particular scaling leads to a very different limit behaviour. We prove that the energy of the effective material depends not only on the strain tensor (as a classical elastic material) but also on the second gradient of the displacement. This result has been announced in [START_REF] Pideri | Un résultat d'homogénéisation pour un matériau élastique renforcé périodiquement par des fibres élastiques de très grande rigidité[END_REF].

Materials whose energy depends on the second gradient of the displacement cannot be considered as Cauchy continua otherwise one would be led to a thermodynamic paradox [START_REF] Dunn | Interstitial working and non classical continuum thermodynamics[END_REF]. This paradox can be removed by extending the thermodynamical framework [START_REF] Dunn | Interstitial working and non classical continuum thermodynamics[END_REF] but the fundamental point is that the Cauchy stress tensor is not sufficient to describe internal forces [START_REF] Seppecher | A numerical study of a moving contact line in Cahn-Hilliard theory[END_REF]. External forces concentrate along any edge of the boundary and the Cauchy theorem defining the Cauchy stress tensor cannot be applied [10,[START_REF] Dell'isola | Edge Contact Forces and Quasi-Balanced Power[END_REF]. Moreover, a supplementary boundary condition is needed to write well-posed problems, which is unusual and not intuitive [START_REF] Seppecher | Etude des conditions aux limites en théorie du second gradient: cas de la capillarité[END_REF]. The simplest way to describe these media is to use the second gradient theory [START_REF] Germain | La méthode des puissances virtuelles en mécanique des milieux continus. Premiere partie: Théorie du second gradient[END_REF][START_REF] Germain | The method of virtual power in continuum mechanics. Part 2: Microstructure[END_REF] or to consider them as Cosserat media [START_REF] Cosserat | Sur la théorie des corps déformables[END_REF]. Our result gives a new example of such a material together with a "microscopic" interpretation of its special features.

We emphasize that, going to the limit, the differential order of the energy changes (as does the system of partial differential equations associated with equilibrium). Such a change is not usual in homogenization theory. It arises in rod or plate theories [START_REF] Acerbi | Thin inclusions in linear elasticity: a variational approach[END_REF] but seems then to be connected with a change of dimension. Our result shows that this is not necessary. Notice also that such a change in the differential order of the energy can not arise when considering scalar problems (like thermal conductivity problems). Indeed, consider a sequence of energies which are quadratic functions of the gradient of a scalar quantity u; these energies decrease when truncating u and this property is preserved when going to the limit. Then, a representation theorem for Dirichlet forms [6, 5] assures that the limit energy can be represented as the sum of a term depending on u and ∇u and a non-local term of the form:

Ω×Ω K (x, y)(u(x ) -u(y)) 2 dxdy.
In other words, we can expect non-local effects but no increase of the differential order. Our result shows that this argument cannot be extended to elasticity problems.

Non-local effects actually arise for some scalar singular perturbation problems [START_REF] Bellieud | Homogénéisation de problèmes elliptiques avec effets non locaux[END_REF][START_REF] Bellieud | Homogénéisation de problèmes elliptiques en présence de fibres de grandes conductivité[END_REF] and we should probably have obtained non-local effects if assuming that ε 2 | log(r ε )| converges to a finite positive value instead of zero. We do not have non-local effects under our assumptions: the second gradient part of the limit energy cannot be interpreted, as it is often done, as the limit of non-local interactions whose range is very short.

Our study is variational. We identify the Γ -limit E 0 of the energy E ε of our composite material. The notion of Γ -convergence corresponds to the intuitive notion of convergence of models: the result is obtained without considering external forces, it remains valid in presence of body forces (for definition and properties of Γ -convergence, refer to [START_REF] Maso | The relationship between edge contact forces and interstitial working allowed by the principle of virtual power[END_REF]).

The limit energy is made explicit in Sect. 2 where we state precisely our result. Section 3 is devoted to the more difficult part of the proof: considering a sequence of displacement fields (u ε ) converging to some u, we have to express the lower bound for the energy E ε (u ε ) in terms of u. This needs an accurate description of the asymptotic behaviour of u ε . Especially u ε has to be described at the scale r ε inside the fibers: we need a multiscale notion of convergence. However, we do not expect any periodicity with period r ε ; the classical notions of multiscale convergence (as defined in [START_REF] Nguetseng | A general convergence result for a functional related to the theory of homogenization[END_REF] or [START_REF] Allaire | homogenization and two scale convergence[END_REF]) are not convenient. In Subsect. 3.1, we develop an adapted notion of double-scale convergence which describes the asymptotic behaviour of u ε in the fibers, that is in a set of scale r ε but with periodicity ε. Section 4 is devoted to the end the proof: for any admissible displacement field u we have to construct an approximating sequence u ε whose limit energy is not larger than E 0 (u). Such an approximation is obtained by choosing u ε = u in the main part of the matrix, a rod-like displacement field in the fibers and a suitable interpolation in transition layers around each fiber.

The main result

Notations and Geometry

In IR 3 we refer to a point x by its Cartesian coordinates (x 1 , x 2 , x 3 ). In the same way the coordinates of any vector u are denoted by (u 1 , u 2 , u 3 ). The symmetric part of the gradient of u (the strain tensor) is denoted by e(u) := (∇u + ∇u t )/2. This tensor belongs to the set of 3-3 symmetric matrices which we denote by M. The trace of a matrix A is denoted by Tr(A).

We use the summation convention, but, as we consider two and three dimensional spaces, we adopt the following convention: a repeated Latin index is summed from 1 to 3 while a Greek index is summed from 1 to 2. In order to describe the periodic structure of our composite material, we introduce two positive real parameters ε and r ε (r ε ≤ ε). Then we define the projection p ε :

p ε : I R 2 -→ IR 2 (x 1 , x 2 ) -→ ε(E( x 1 ε ) + 1 2 ), ε(E( x 2 ε ) + 1 2 )
,

where E(t) denotes the integer part of a real t and we define the periodic function y ε by

y ε : I R 3 -→ IR 2 (x 1 , x 2 , x 3 ) -→ r -1 ε [(x 1 , x 2 ) -p ε (x 1 , x 2 )] .
Next, we define the sets F ε and M ε , referred to as "the fibers" and "the matrix" respectively, by:

F ε := {x ∈ IR 3 : |y ε (x )| < 1} , M ε := IR 3 \ F ε (1)
We assume that the composite material lies in the cube Ω = (0, 1) 3 and we denote by B its "lower" face: B = (0, 1) 2 × {0}. On Ω the projection p ε ranges onto a finite set of points which we denote by x p ; p ∈ P ε := {1 . . . 1 ε 2 } . The domain Ω is the union of the ε -2 parallelepipeds P p ε := {x ∈ Ω : p ε (x ) = x p } which correspond to the periods of y ε .

The fiber contained in the period P p ε (a circular cylinder of radius r ε , see Fig. 2) is denoted by

F p ε := F ε ∩P p ε . The total volume of the fibers contained in Ω is |F ε ∩ Ω| = p∈∩P ε |F p ε | = πr 2 ε ε -2 .

Elastic energy

We assume that F ε and M ε consist of two different isotropic elastic materials: we define for every Borel set D the matrix energy by

E m (D, u) := D [ λ 0 2 (Tr(e(u))) 2 + µ 0 e(u) 2 ] dx , (2)
where (λ 0 , µ 0 ) denote the (positive) elasticity Lamé coefficients in the matrix. In the same way, we define the fiber energy by

E f ε (D, u) := D [ λ ε 2 (Tr(e(u))) 2 + µ ε e(u) 2 ] dx , (3)
where (λ ε , µ ε ) denote the (positive) elasticity Lamé coefficients in the fibers.

We assume perfect adhesion between the matrix and the fibers. Moreover, we assume that both materials are fixed to the plane {x 3 = 0}. Then, for any displacement u ∈ L 2 (Ω, IR 3 ), we define the total energy 3 ) and u = 0 on B , +∞ otherwise .

E ε (u) := E ε (Ω, u) of our composite material by E ε (Ω, u) :=    E m (M ε ∩ Ω, u) + E f ε (F ε ∩ Ω, u), if u ∈ H 1 (Ω, IR
(4)

Order of magnitude of the different parameters

In order to study the Γ -limit of E ε as ε tends to zero, we must specify our assumptions upon the behaviour of r ε , λ ε and µ ε as ε tends to zero: we assume that r ε obeys the limit relations

lim ε-→0 r ε ε = 0 , (5)
lim

ε-→0 ε 2 log(r ε ) = 0 (6)
and that µ ε and λ ε fulfill the limit conditions

lim ε-→0 µ ε r 4 ε ε 2 = µ 1 ∈]0, +∞[ , (7) lim
ε-→0 λ ε µ ε = ∈ [0, +∞[ . (8)
Assumption [START_REF] Bellieud | Homogénéisation de problèmes elliptiques avec effets non locaux[END_REF] states that the fibers are much thinner than the period of the medium; it is one of our basic assumptions. However, they cannot be too thin, otherwise the connection between the displacement fields in the matrix and in the fibers disappears when ε tends to zero. This fact can be explained as follows: if the radius of the fibers is infinitely smaller than ε, the fibers behave like one-dimensional media and it is well known that a one-dimensional medium has no connection with a three-dimensional elastic medium. We will see later that restriction (6) assures that the global displacement of each fiber coincides with the displacement of the matrix when ε tends to zero. Note that assumption (6) is not very restrictive: any power law r ε = ε γ , γ > 1 is admissible.

The energy of a bent rod is related to its curvature, that is to the second gradient of its displacement. As our goal is to obtain second gradient effects, we expect each fiber to behave like a rod. The bending stiffness of a unique fiber is

π 4 r 4 ε µ ε 3λ ε + 2µ ε λ ε + µ ε
(refer to any textbook for mechanics of structures or to the pioneering work of St. Venant [START_REF] De Saint-Venant | Mémoire sur la flexion des prismes[END_REF]). Assumptions ( 7) and ( 8) state that this stiffness is of the order of ε 2 , the inverse of the number (ε -2 ) of fibers.

The main result

Our result states that E ε Γ -converges in L 2 (Ω, IR 3 ) to E 0 defined by:

E 0 (u) =                    E m (Ω, u) + Ω k 2 ( ∂ 2 u 1 ∂x 2 3 ) 2 + ( ∂ 2 u 2 ∂x 2 3 ) 2 dx , if u ∈ H 1 (Ω, IR 3 ), ∂ 2 u ∂x 2 3 ∈ L 2 (Ω, IR 3 ), u 3 = 0 a.e. in Ω, u = ∂u ∂x 3 = 0 a.e. on B , +∞ otherwise. ( 9
)
where

k = π 4 3 + 2 + 1 µ 1 , (10)
More precisely we have the following:

Theorem 1. i) Let u ε be a sequence such that E ε (u ε ) is bounded. Then u ε is strongly relatively compact in L 2 (Ω, IR 3 ).
ii) Moreover, for any sequence u ε converging to u in L 2 (Ω, IR 3 ), the following lower bound inequality holds:

lim inf ε-→0 E ε (u ε ) ≥ E 0 (u). ( 11
)
iii) Conversely, for every u in L 2 (Ω, IR 3 ), there exists an approximating sequence u ε in L 2 (Ω, IR 3 ) such that

u ε -→ u in L 2 (Ω, IR 3 ), lim sup ε-→0 E ε (u ε ) ≤ E 0 (u). ( 12
)
Proof of assertion (i): It is clear from assumptions (5), ( 7) and ( 8) that λ ε and µ ε tend to infinity. Then there exists a positive real c such that E ε (u) ≥ c Ω e(u) 2 dx for every u in H 1 (Ω, IR 3 ). Due to Korn's inequality, there exists a positive real C such that E ε (u) ≥ C ||u|| H 1 (Ω,IR 3 ) . The sequence u ε is then bounded in H 1 (Ω, IR 3 ): it is strongly relatively compact in L 2 (Ω, IR 3 ).

The proofs of (ii) and (iii) are less straightforward. They are given in the following two sections.

Proof of the lower bound inequality

Preliminaries, double-scale convergence

Let us denote by D 1 the unit disk of IR 2 and by D the set of functions D := C ∞ c (Ω × D 1 , IR). We associate to the sequence of sets (F ε ∩ Ω) the following "double scale" convergence:

Definition: We say that a sequence u ε in L 2 (Ω, IR) double scale converges to v ∈ L 2 (Ω × D 1 , IR) and we write u ε v if and only if ∀ϕ ∈ D , - Fε∩Ω u ε (x ) ϕ(x , y ε (x )) dx → - Ω - D1 v(x, y) ϕ(x, y) dydx . ( 13
)
This definition is extended to vector field or tensor field sequences: we say that such sequences d.s.-converge if and only if every component is d.s.-convergent.

Remark 1. For every function

Φ ∈ D , u ε v =⇒ Φ(., y ε (.)) u ε (.) Φ v . (14)
Indeed, for every ψ ∈ D , the product ψ Φ belongs to D and the result is obtained by applying the definition of the d.s.-convergence of u ε with ϕ = ψΦ.

Lemma 1. For every function Φ ∈ D we have

Φ(., y ε (.)) Φ . ( 15
)
Proof: For ϕ ∈ D , let us compute the limit of -Fε∩Ω ϕ(x , y ε (x ))dx . Using the Fubini theorem and changing variables in each fiber we get

- Fε∩Ω ϕ(x , y ε (x ))dx = ε 2 p∈Pε - F p ε ϕ(x , y ε (x )) dx = |F ε ∩ Ω| -1 r 2 ε p∈Pε 1 0 D1 ϕ((x p 1 + r ε y 1 , x p 2 + r ε y 2 , x 3 ), (y 1 , y 2 )) dy dx 3 .
As the function ϕ is uniformly continuous on Ω × D 1 , we have the following uniform estimations:

|ϕ((x p 1 + r ε y 1 , x p 2 + r ε y 2 , x 3 ), (y 1 , y 2 )) -ϕ((x p 1 , x p 2 , x 3 ), (y 1 , y 2 ))| = O(r ε ) and p∈Pε ϕ((x p 1 , x p 2 , x 3 ), y)1 P p ε (x ) -ϕ(x , y) = O(ε) , which implies 1 0 ε 2 p∈Pε ϕ((x p 1 , x p 2 , x 3 ), y) dx 3 - Ω ϕ(x , y) dx = O(ε). Hence lim ε→0 - Fε∩Ω ϕ(x , y ε (x )) dx = lim ε→0 |F ε ∩ Ω| -1 r 2 ε ε -2 Ω D1 ϕ(x , y) dy dx = - Ω - D1 ϕ(x , y) dy dx .
In other words, the constant function 1 d.s.-converges to itself. The lemma is proved by recalling Remark 1.

Lemma 2. Let u ε be a sequence in L 2 (Ω, IR) such that -Fε∩Ω u 2 ε (x ) dx is bounded, then there exists a subsequence of u ε (still denoted by u ε ) and a function v

∈ L 2 (Ω × D 1 , IR) such that u ε v Proof: Assume -Fε∩Ω u 2 ε (x ) dx ≤ M
and consider the sequence of measures ν ε on Ω × D 1 defined by

ν ε := |F ε ∩ Ω| -1 u ε (x ) δ yε(x ) (dy) dx . (16)
Since the sequence ν ε (Ω × D 1 ) is bounded, there exists a measure ν such that ν ε ν for some subsequence. Moreover, for every ϕ ∈ D , we have

Ω×D1 ϕ(x , y)d ν ε = |F ε ∩ Ω| -1 Fε∩Ω ϕ(x , y ε (x ))u ε (x )dx ≤ - Fε∩Ω (u ε (x )) 2 dx 1/2 - Fε∩Ω (ϕ(x , y ε (x ))) 2 dx 1/2 ≤ M 1/2 - Fε∩Ω ϕ(x , y ε (x )) 2 dx 1/2 .
As ϕ 2 ∈ D , using Lemma 1, we have lim sup

ε→0 Ω×D1 ϕ(x , y)d ν ε ≤ M 1/2 - Ω - D1 (ϕ(x , y)) 2 dydx 1/2 Ω×D1 ϕ(x , y)d ν ≤ M 1/2 π -1/2 ||ϕ|| L 2 (Ω×D1) .
The measure ν, as a linear functional, is bounded on the unit ball of L 2 (Ω × D 1 , IR): there exists a function v ∈ L 2 (Ω × D 1 , IR) such that ν = v dx dy. The convergence of the sequence of measures ν ε to the measure v dx dy is clearly equivalent to the d.s.-convergence of u ε to v.

Let us notice that Lemma 2 can obviously be extended to vector or matrix fields.

Lemma 3. Let u ε be a bounded sequence in H 1 (Ω, IR 3 ). Then, there exists a constant C such that, for ε small enough,

- Fε∩Ω (u ε (x )) 2 dx ≤ C ||u ε || L 2 (Ω,IR 3 ) -ε 2 log(r ε ) . (17)
Proof: Assume that ||u ε || 2 H 1 (Ω,IR 3 ) ≤ M . Then

p∈Pε P p ε (∇u ε ) 2 dx ≤ M .
In each period P p ε , we use the cylindrical coordinates, defining u p ε by

u p ε (r, θ, x 3 ) := u ε (x p 1 + r cos θ, x p 2 + r sin θ, x 3 ) . (18)
Then, we have, for every

ρ 1 ≤ ρ 2 ≤ ε/2, p∈Pε 1 0 2π 0 ρ2 ρ1 ∂u p ε ∂r 2 rdrd θdx 3 ≤ M .
A simple one-dimensional minimization shows that ρ2 ρ1

( ∂u p ε ∂r ) 2 rdr ≥ [u p ε (ρ 2 ) -u p ε (ρ 1 )] 2 log(ρ 2 ) -log(ρ 1 )
.

Hence

p∈Pε 1 0 2π 0 [u p ε (ρ 2 , θ, x 3 ) -u p ε (ρ 1 , θ, x 3 )] 2 d θdx 3 ≤ M log( ρ 2 ρ 1
).

Let us denote by f the quantity

f (ρ) = p∈Pε 1 0 2π 0 (u p ε (ρ, θ, x 3 )) 2 d θ dx 3 . (19)
The last inequality implies that, for every

ρ 1 ≤ ρ 2 ≤ ε/2, f (ρ 1 ) ≤ 2f (ρ 2 ) + 2M log( ρ 2 ρ 1
) .

As the ratio r ε /ε tends to zero, we may assume, without loss of generality, that r ε ≤ ε/4. Then, for every

ρ 2 ∈ [ε/4, ε/2], we can bound the mean value of u 2 ε on F ε ∩ Ω by - Fε∩Ω u 2 ε dx = |F ε ∩ Ω| -1 rε 0 f (r) r dr ≤ |F ε ∩ Ω| -1 rε 0 2r f (ρ 2 ) + 2Mr log( ρ 2 r ) dr ≤ ε 2 πr 2 ε f (ρ 2 ) r 2 ε + Mr 2 ε (log( ρ 2 r ε ) + 1 2 ) ≤ 1 π 4ε ρ 2 f (ρ 2 ) + M ε 2 (log( ε r ε ) + 1 2 )
and, taking the mean value of this last term for

ρ 2 ∈ [ε/4, ε/2], we get - Fε∩Ω u 2 ε dx ≤ 1 π 16 ε 2 ε 4 f (r) r dr + M ε 2 (log(ε) -log(r ε ) + 1 2 ) ≤ 1 π 16||u ε || 2 L 2 (Ω,IR 3 ) + M ε 2 (log(ε) -log(r ε ) + 1 2 )
.

For ε sufficiently small, | log(ε) + 1/2| ≤ | log(r ε )|. The lemma is proved by taking C = sup{16/π, 2 M /π}. Lemma 4. Let u ε be a bounded sequence in H 1 (Ω, IR 3 ). Then i) - Fε∩Ω u 2 ε dx is bounded. ii) If u ε → u in L 2 (Ω, IR 3 ), then - Fε∩Ω (u ε (x ) -u(x )) 2 dx → 0. iii) If u ε → u in L 2 (Ω, IR 3 ) and u ε v, then u(x ) = - D1 v(x, y) dy , a.e. in Ω. (20)
Proof: Assertion (i) is a trivial consequence of Lemma 3 and assumption (6). Here it becomes clear how assumption (6) connects the displacement in the fibers to the displacement in the matrix. Note that, at this point, the boundedness of ε 2 log(r ε ) should be sufficient. Assertion (ii) needs the convergence ε 2 log(r ε ) → 0. Then one simply must apply Lemma 3 to the sequence (u εu).

To prove assertion (iii), let us consider for any ν > 0, a field

Φ ν ∈ C ∞ c (Ω, IR 3 ) such that ||Φ ν - u|| L 2 (Ω,IR 3 ) < ν. For any ϕ ∈ C ∞ c (Ω, IR 3 ) we have | lim ε→0 - Fε∩Ω (u ε (x ) -Φ ν (x ))ϕ(x )dx | ≤ lim ε→0 - Fε∩Ω |u ε (x ) -Φ ν (x )| 2 dx 1/2 lim ε→0 - Fε∩Ω ϕ(x ) 2 dx 1/2
.

Applying Lemma 3 to the sequence |u ε -Φ ν | shows that this last term is bounded by the norm ||u ε (x ) -Φ ν (x )|| L 2 (Ω,IR 3 ) and therefore is of order O(ν). Now, passing to the double scale limit, using the definition of v and Lemma 1, we get

| - Ω (- D1 v(x, y) dy -Φ ν (x )) ϕ(x ) dx | ≤ O(ν) .
Assertion (iii) is proved by recalling that this inequality is valid for every ν.

Limits of a sequence with bounded energy

Lemma 5. Let u ε be a sequence of L 2 (Ω, IR 3 ) with bounded energy. Then, up to a subsequence (still denoted by u ε ), there exist

v ∈ L 2 (Ω × D 1 , IR 3 ), w ∈ L 2 (Ω × D 1 , IR) and χ ∈ L 2 (Ω × D 1 , M) such that u ε v , u ε3 r ε w , e(u ε ) r ε χ . ( 21
)
Proof: Assume E ε (u ε ) < M , then the sequence u ε is bounded in H 1 (Ω, IR 3 ), Lemma 4 states that the sequence

- Fε∩Ω u 2 ε dx is bounded and Lemma 2 implies the existence of v ∈ L 2 (Ω × D 1 , IR 3 ) such that u ε v.
On the other hand, as µ ε Fε∩Ω e(u ε ) 2 dx < M , we have

µ ε r 4 ε ε 2 π - Fε∩Ω ( 1 r ε ∂u ε 3 ∂x 3 ) 2 dx < M .
As any sequence with bounded energy satisfies u ε 3 (x 1 , x 2 , 0) = 0 a.e. on B , a simple one-dimensional minimization shows that 1 0

( ∂u ε 3 ∂x 3 ) 2 dx 3 ≥ π 2 4 1 0 (u ε 3 ) 2 dx 3 , for a.e. (x 1 , x 2 ).
Hence, 2 dx is bounded: the sequence u ε 3 /r ε satisfies the assumptions of Lemma 2; the existence of w is assured.

µ ε r 4 ε ε 2 - Fε∩Ω ( u ε 3 r ε ) 2 dx ≤ 4 π 3 M . As µ ε r 4 ε /ε 2 → µ 1 , the sequence - Fε∩Ω ( u ε 3 r ε )
In the same way, from inequality µ ε Fε∩Ω e ( u ε ) 2 dx < M , we deduce

µ ε r 4 ε ε 2 π - Fε∩Ω ( e(u ε ) r ε ) 2 dx < M .
The sequence e(u ε )/r ε verifies the assumptions of Lemma 2: the existence of χ is assured.

Lemma 6. Consider a sequence u ε with bounded energy and converging to some u in L 2 (Ω, IR 3 ), then

u ∈ H 1 (Ω, IR 3 ) , ∂ 2 u 1 ∂x 2 3 ∈ L 2 (Ω, IR 3 ) , ∂ 2 u 2 ∂x 2 3 ∈ L 2 (Ω, IR 3 ) , u 3 (x ) = 0 a.e. in Ω .
Moreover, there exists a subsequence (still denoted by u ε ) and q ∈ L 2 (Ω, IR) such that e(u ε )

r ε 33 q(x ) - ∂ 2 u 1 ∂x 2 3 (x ) y 1 - ∂ 2 u 2 ∂x 2 3 (x ) y 2 . ( 22
)
Proof: First, let us notice that the sequence u ε is bounded in H 1 (Ω, IR 3 ). Then the limit u belongs to H 1 (Ω, IR 3 ). Lemma 5 assures the existence of v ∈ L 2 (Ω × D 1 , IR 3 ), w ∈ L 2 (Ω × D 1 , IR) and χ ∈ L 2 (Ω × D 1 , M) such that, up to a subsequence,

u ε v , u ε3 r ε w and e(u ε ) r ε χ .
The convergence u ε3 /r ε w immediately yields u ε3 0, i.e. v 3 = 0. Using the relation u(x ) = -D1 v(x, y) dy stated in Lemma 4, we get the identity u 3 = 0 a.e. in Ω.

Consider now a tensor field ϕ ∈ C ∞ c (Ω × D 1 , M). We have, using the definition of χ and the divergence theorem,

- Ω - D1 χ ij (x , y) ϕ ij (x , y) dy dx = lim ε→0 - Fε∩Ω 1 r ε e ij (u ε )(x ) ϕ ij (x , y ε (x )) dx , = lim ε→0 - Fε∩Ω 1 r ε ∂u ε i ∂x j (x ) ϕ ij (x , y ε (x )) dx , = -lim ε→0 - Fε∩Ω 1 r ε u ε i (x ) ∂ϕ ij ∂x j (x , y ε (x )) + 1 r ε ∂ϕ i α ∂y α (x , y ε (x )) dx , = -lim ε→0 1 r 2 ε - Fε∩Ω u ε β (x) ∂ϕ βα ∂y α (x , y ε (x )) dx + 1 r ε - Fε∩Ω u ε β (x) ∂ϕ βj ∂x j (x , y ε (x )) + u ε 3 (x ) r ε ∂ϕ 3α ∂y α (x , y ε (x )) dx + - Fε∩Ω u ε 3 (x ) r ε ∂ϕ 3j ∂x j (x , y ε (x )) dx . ( 23
)
Multiplying equation ( 23) by r 2 ε and passing to the limit ε → 0 gives

lim ε→0 - Fε∩Ω u ε β (x) ∂ϕ βα ∂y α (x , y ε (x )) dx = 0 , - Ω - D1 v β (x, y) ∂ϕ βα ∂y α (x , y) dy dx = 0 , ∂v β ∂y α ϕ βα = 0 , ( 24
)
where <> denotes the distribution bracket on Ω × D 1 . This last equation, valid for any field ϕ of a symmetric plane matrix and whose support is included in Ω × D 1 , is equivalent to the antisymmetry (in the sense of distributions)

∂v 1 ∂y 2 = - ∂v 2 ∂y 1 , ∂v 1 ∂y 1 = ∂v 2 ∂y 2 = 0 .
Then (refer for instance to [START_REF] Dret | Problèmes variationnels dans les multi-domaines[END_REF]) there exist three functions c 1 , c 2 and t in L 2 (Ω, IR) such that

v 1 (x , y) = c 1 (x) -t(x) y 2 , v 2 (x,y) = c 2 (x) + t(x) y 1 .
Lemma 4 implies c 1 = u 1 and c 2 = u 2 . Hence

v 1 (x , y) = u 1 (x) -t(x) y 2 , v 2 (x,y) = u 2 (x) + t(x) y 1 . ( 25
)
Now, consider the fields ϕ such that ϕ βα = 0, ∀ α, β ∈ {1, 2}. Multiplying equation ( 23) by r ε and passing to the limit gives The Schwarz theorem implies that ∂t/∂x 3 = 0; then

lim ε→0 - Fε∩Ω u ε β (x) ∂ϕ β3 ∂x 3 (x , y ε (x )) + u ε 3 (x ) r ε ∂ϕ 3α ∂y α (x , y ε (x )) dx = 0 , - Ω - D1 v β (x,
- ∂u 1 ∂x 3 = ∂w ∂y 1 , - ∂u 2 ∂x 3 = ∂w ∂y 2 .
Therefore there exists a function s in L 2 (Ω, IR) such that

w(x , y) = - ∂u α ∂x 3 y α + s(x) . (27)
Finally, considering matrix fields ϕ with a unique non vanishing component ϕ 33 , equation ( 23) leads to

- Ω - D1 χ 33 (x , y) ϕ 33 (x , y) dx = -lim ε→0 - Fε∩Ω u ε 3 r ε (x ) ∂ϕ 33 ∂x 3 (x , y ε (x )) dx = -- Ω - D1
w(x, y) ∂ϕ 33 ∂x 3 (x , y) dy dx

χ 33 ϕ 33 = ∂w ∂x 3 ϕ 33 . ( 28
)
Then χ 33 = ∂w/∂x 3 in the sense of distributions. As χ belongs to L 2 (Ω × D 1 , M), ∂w/∂x 3 belongs to L 2 (Ω × D 1 , IR). This means, by using (27) that ∂ 2 u α /∂x 2 3 ∈ L 2 (Ω, IR), q := ∂s/∂x 3 ∈ L 2 (Ω, IR) and

χ 33 (x , y) = - ∂ 2 u α ∂x 2 3 (x) y α + q(x) . ( 29
)

Lower bound for the energy

Let u ε be a sequence with bounded energy converging to some u in L 2 (Ω, IR 3 ). We can assume without loss of generality that E ε (u ε ) converges to lim inf E ε (u ε ). Then assertion (ii) of Theorem 1 will be proved if we prove that for some subsequence (still denoted by u ε ) we have

lim inf ε→0 E ε (u ε ) ≥ E 0 (u) .
First, let us recall that the sequence u ε is bounded in H 1 (Ω, IR 3 ), then

u ∈ H 1 (Ω, IR 3 ) . (30)
It is easy to get the lower bound for the energy outside the fibers: indeed, as

E ε (u ε ) is bounded, E f ε (F ε ∩ Ω, u ε ) is also bounded. As the ratios µ 0 /µ ε and λ 0 /λ ε tend to zero, then E m (F ε ∩ Ω, u ε ) tends to 0. Hence lim inf ε→0 E m (M ε ∩ Ω, u ε ) = lim inf ε→0 E m (Ω, u ε ) ≥ E m (Ω, u) . (31)
To estimate the energy in the fibers we use the lemmas stated in the preceeding subsections. Indeed we have

lim inf ε→0 E f ε (F ε ∩ Ω, u ε ) = lim inf ε→0 µ ε r 4 ε ε 2 π - Fε∩Ω ( e(u ε ) r ε ) 2 + λ ε µ ε ( Tr(e(u ε )) r ε ) 2 dx ≥ πµ 1 lim inf ε→0 - Fε∩Ω ( e(u ε ) r ε ) 2 + 2 ( Tr(e(u ε )) r ε ) 2 dx . (32)
From Lemma 5, we know that, possibly passing to a subsequence, e(u ε )/r ε admits a double scale limit χ.

As we cannot pass to the limit directly in inequality (32), we write its dual form

lim inf ε→0 E f ε (F ε ∩ Ω, u ε ) ≥ sup ϕ πµ 1 lim inf ε→0 - Fε∩Ω e(u ε (x )) r ε : ϕ(x , y ε (x ))- 1 4 ϕ(x , y ε (x )) 2 + 4(2 + 3 ) (Tr(ϕ(x , y ε (x )))) 2 dx ,
where the supremum is taken for every ϕ ∈ C ∞ c (Ω × D 1 , M). Then Remark 1 and Lemma 1 allow to pass to the limit

lim inf ε→0 E f ε (F ε ∩ Ω, u ε ) ≥ sup ϕ πµ 1 - Ω - D1 χ(x , y) : ϕ(x, y)- 1 4 (ϕ(x, y)) 2 + 4(2 + 3 ) (Tr(ϕ(x
, y))) 2 dy dx .

As C ∞ c (Ω × D 1 , M) is dense in L 2 (Ω × D 1 , M), we get lim inf ε→0 E f ε (F ε ∩ Ω, u ε ) ≥ πµ 1 - Ω - D1 χ 2 (x , y) + 2 (Tr(χ(x , y))) 2 dy dx .
It is easy to verify that, for every M in M,

M 2 + 2 (Tr(M )) 2 ≥ 3 + 2 2( + 1) M 2 33 . (33)
Hence,

lim inf ε→0 E f ε (F ε ∩ Ω, u ε ) ≥ πµ 1 3 + 2 ( + 1) - Ω - D1 χ 2 33 (x , y)dy dx . (34)
From Lemma 6, we know that

∂ 2 u α ∂x 2 3 ∈ L 2 (Ω, IR) , (35)
and we can express χ 33 in terms of these second derivatives of u Then we may deduce that

lim inf ε→0 E f ε (F ε ∩ Ω, u ε ) ≥ π 2 µ 1 3 + 2 ( + 1) - Ω - D1 q(x) - ∂ 2 u α ∂x 2
lim inf ε→0 E f ε (F ε ∩ Ω, u ε ) ≥ πµ 1 (3 + 2) 8 ( + 1) - Ω 4 q 2 (x) + ∂ 2 u 1 ∂x 2 3 2 + ∂ 2 u 2 ∂x 2 3 2 dx , which implies lim inf ε→0 E f ε (F ε ∩ Ω, u ε ) ≥ k 2 - Ω ( ∂ 2 u 1 ∂x 2 3 ) 2 + ( ∂ 2 u 2 ∂x 2 3 ) 2 dx , (36)
where k is defined by (10).

In order to obtain the boundary conditions, let us consider the extended domain Ω := (0, 1) 2 ×] -1, 1[ and the extensions ũε and ũ of u ε and u on Ω defined by ũε := u ε on Ω , ũ := u on Ω , ũε := ũ := 0 on Ω \ Ω .

The sequence E ε ( Ω, ũε ) is bounded and ũε converges to ũ in L 2 ( Ω, IR 3 ); thus the results of Lemma 6 can be applied: ũ ∈ H 1 ( Ω, IR 3 ) and ∂ 2 ũα /∂x || L 2 (Ω,IR 3 ) .

For any u ∈ L 2 (Ω, IR 3 ) such that E 0 (u) < +∞, i.e., for any u ∈ H , we have to construct an approximating sequence u ε in L 2 (Ω, IR 3 ) such that u ε -→ u in L 2 (Ω, IR 3 ) and lim sup

ε-→0 E ε (u ε ) ≤ E 0 (u).
It is easy to verify that

H := {u ∈ C ∞ (Ω, IR 3 ), u = ∂u ∂x 3 = ∂ 2 u ∂x 2 3 = 0 a.e
C ε := {x ∈ Ω : 1 < |y ε (x)| < r -1 ε R ε }, B ε := {x ∈ Ω : |y ε (x )| > r -1 ε R ε },
The part of C ε contained in a period P p ε is denoted by C p ε := C ε ∩ P p ε For every p in P ε , we define the function v p ε ∈ C ∞ ((0, 1), IR 3 ) by

v p ε (x 3 ) := - D1 u(x p 1 + r ε y 1 , x p 2 + r ε y 2 , x 3 ) dy 1 dy 2 (38)
and the function w p ε ∈ C ∞ ((0, 1) × IR 2 , IR 3 ) by

w p ε1 (x 3 , y) := v p ε1 (x 3 ) + r 2 ε 2( + 1) ∂ 2 v p ε1 ∂x 2 3 y 2 1 -y 2 2 2 + ∂ 2 v p ε2 ∂x 2 3 y 1 y 2 , w p ε2 (x 3 , y) := v p ε2 (x 3 ) + r 2 ε 2( + 1) ∂ 2 v p ε2 ∂x 2 3 y 2 2 -y 2 1 2 + ∂ 2 v p ε1 ∂x 2 3 y 1 y 2 , (39)
w p ε3 (x 3 , y) := -r ε ∂v p εα ∂x 3 y α .
The function w p ε may be interpreted as the rod-like displacement of the fiber F p ε whose global displacement is v p ε [START_REF] De Saint-Venant | Mémoire sur la flexion des prismes[END_REF]. As u ∈ H , we have u = ∂u/∂x 3 = ∂ 2 u/∂x 2 3 = 0 on B . Therefore every fonction w p ε vanishes for x 3 = 0.

We define now the approximating sequence (u ε ) by setting

u ε (x ) :=              u(x) o n B ε , w p ε ( x 3 , y ε ( x )) on each fiber F p ε , γ(r) w p ε (x 3 , (cos θ, sin θ)) on each transition +(1 -γ(r)) u(x ) layer C p ε , (40)
where (r, θ) denote the polar coordinates defined in each period P p ε by x 1 = x p 1 + r cos θ, x 2 = x p 2 + r sin θ and γ is the function defined by

γ(r) := log(r) -log(R ε ) log(r ε ) -log(R ε ) .
Notice that, by construction, u ε belongs to H 1 (Ω, IR 3 ) and satisfies u ε = 0 on B . Then

E ε (u ε ) = E m (B ε ∩ Ω, u ε ) + E m (C ε ∩ Ω, u ε ) + E f ε (F ε ∩ Ω, u ε ) . (41)
Moreover, u ε tends to u in L 2 (Ω, IR 3 ): indeed u ε coincides with u on B ε , |Ω \ B ε | → 0, and (u ε ) is uniformly bounded on F ε and C ε .

Estimation for the energy of u ε in the matrix

As u ε (x ) := u(x) on B ε , we have

E m (B ε , u ε ) = E m (B ε , u). As R ε /ε → 0 one has |Ω \ B ε | → 0. Moreover, u ∈ H 1 (Ω), then E m (Ω \ B ε , u) → 0 and lim ε→0 E m (B ε , u ε ) = E m (Ω, u) . ( 42
)

Estimation for the energy of u ε in the fibers

Let us estimate the energy of u ε in each fiber

F p ε : As u ε (x) = w p ε (x, y ε (x)) in F p ε , we have e 11 (u ε ) = e 22 (u ε ) = r ε 2( + 1) ∂ 2 v p εα ∂x 2 3 y α , e 33 (u ε ) = -r ε ∂ 2 v p εα ∂x 2 3 y α , e 12 (u ε ) = e 21 (u ε ) = 0 , e 13 (u ε ) = e 31 (u ε ) = r 2 ε 4( + 1) ∂ 3 v p ε1 ∂x 3 3 y 2 1 -y 2 2 2 + ∂ 3 v p ε2 ∂x 3 3 y 1 y 2 , e 23 (u ε ) = e 32 (u ε ) = r 2 ε 4( + 1) ∂ 3 v p ε2 ∂x 3 3 y 2 2 -y 2 1 2 + ∂ 3 v p ε1 ∂x 3 3 y 1 y 2 . Hence E f ε (F p ε , u ε ) = r 4 ε 1 0 D1 µ ε 3 2 + 4 + 2 2( + 1) 2 + λ ε 1 2( + 1) 2 ∂ 2 v p εα ∂x 2 3 y α 2 dx +r 6 ε 1 0 D1 µ ε 2 8( + 1) 2 ∂ 3 v p ε1 ∂x 3 3 2 + ∂ 3 v p ε2 ∂x 3 3 2 y 2 1 + y 2 2 2 2 dx .
Computing the integrals on D 1 and summing for all sets F p ε we get

E f ε (F ε ∩ Ω, u ε ) = p∈Pε 2 πr 6 ε µ ε 96( + 1) 2 1 0 ∂ 3 v p ε1 ∂x 3 3 2 + ∂ 3 v p ε2 ∂x 3 3 2 dx 3 + 3 2 + 4 + 2 + λε µε 2( + 1) 2 πr 4 ε µ ε 4 1 0 ∂ 2 v p ε1 ∂x 2 3 2 + ∂ 2 v p ε2 ∂x 2 3 2 dx 3 .
Passing to the limit ε → 0, we have

ε -2 2 8( + 1) 2 πr 6 ε µ ε 12 → 0 , ε -2 3 2 + 4 + 2 + λε µε 2( + 1) 2 πr 4 ε µ ε 4 → k 2 ,
where k is defined by (10). Moreover, using the definition of the functions v p ε , we have

ε 2 p∈Pε 1 0 ∂ 2 v p εα ∂x 2 3 2 dx 3 = ε 2 - Fε∩Ω ∂ 2 u α ∂x 2 3 2 dx , ε 2 p∈Pε 1 0 ∂ 3 v p εα ∂x 3 3 2 dx 3 = ε 2 - ∂ 3 u α ∂x 3 3 2 dx . Hence lim ε→0 E f ε (F ε ∩ Ω, u ε ) = k 2 - Ω ∂ 2 u 1 ∂x 2 3 2 + ∂ 2 u 2 ∂x 2 3 2 dx . (43)

Estimation for the energy of u ε in the transition layer

Let M = sup Ω sup(∇u, ∇ 2 u, ∇ 3 u) . We restrict attention to a cylinder C p ε and prove, in a first step, that |∇u ε | is bounded on C p ε . We use the cylindrical coordinates (r, θ, x 3 ) defined by

x 1 = x p 1 + r cosθ, x 2 = x p 2 + r sinθ (on C p ε we have r ∈ [r ε , R ε ]).
Clearly, in view of the definition of w p ε , there exists a positive real M 1 such that

| ∂w p ε ∂θ (x 3 , (cos θ, sin θ))| ≤ M 1 r ε . Moreover, | ∂u ∂θ | ≤ Mr and whence | 1 r ∂u ε ∂θ | = 1 r |(1 -γ(r)) ∂u ∂θ (x p 1 + r cos θ, x p 2 + r sin θ, x 3 ) + γ(r) ∂w p ε ∂θ (x 3 , (cos θ, sin θ))| ≤ M 1 + M . (44)
On the other hand, owing to the definition of w p ε , there exists a positive real M 2 such that, for every y ∈ D 1 and

x 3 ∈ [0, 1], |w p ε (x 3 , y) -v p ε (x 3 )| ≤ M 2 r ε . From the definition of v p ε , we have, for every r ≥ r ε |v p ε (x 3 ) -u(x p 1 + r cos θ, x p 2 + r sin θ, x 3 )| ≤ 2 M r .
Then, there exists a positive real M 3 such that, for every r ≥ r ε , y ∈ D 1 , θ ∈ [0, 2π] and x 3 ∈ [0, 1],

|w p ε (x 3 , y) -u(x p 1 + r cos θ, x p 2 + r sin θ, x 3 )| ≤ M 3 r .
(C ε , u ε ) = 0 . (47)
Assertion (iii) of Theorem 1 is proved by the estimations (42), ( 43) and (47).

Comments

Due to the properties of Γ -convergence, our result is still valid when external body forces are present. Indeed, a term Ω f (x )u(x )dx can be added to both E ε and E 0 . In that way, we can solve non-trivial equilibrium problems.

Our result states that the homogenized material is a second gradient material: it has a "three dimensional bending stiffness" k . This is not so surprising: it is well know that elastic cylinders, when their radius tends to zero, behave like rods (which are second gradient one-dimensional media): in a sense, we studied the homogenized properties of a system of rods connected by an elastic matrix. However, it must be emphasized that such a result could not be reached by considering directly an elastic matrix reinforced by one-dimensional rods (there is no interaction between a one-dimensional and an elastic three-dimensional medium).

The limit energy E 0 contains a remaining classic elastic part, E m (Ω, u). One could consider, afterwards, the limit (µ 0 , λ 0 ) → (0, 0) in E 0 and obtain an energy depending only on the second gradient of the displacement (the bending stiffness k does not depend on µ 0 or λ 0 ).

The particular features of second gradient materials, like the hyperstress tensor [START_REF] Germain | La méthode des puissances virtuelles en mécanique des milieux continus. Premiere partie: Théorie du second gradient[END_REF][START_REF] Germain | The method of virtual power in continuum mechanics. Part 2: Microstructure[END_REF], flux of interstitial working [START_REF] Dunn | Interstitial working and non classical continuum thermodynamics[END_REF]10], edge forces [START_REF] Dell'isola | Edge Contact Forces and Quasi-Balanced Power[END_REF], presence of a force distribution of order one with respect to the normal derivative [START_REF] Seppecher | Etude des conditions aux limites en théorie du second gradient: cas de la capillarité[END_REF] can be interpreted in our particular case as limits of some microscopic elastic forces.

An open question raised by our study is the general condition for the change of differential order of the energy when passing to the limit. We already pointed out that such a change was impossible for scalar problems. Our feeling is that the properties of the kernel of the energy density (rigid motions in our case) is essential: it leads to constraints verified by the limit of sequences with bounded energy (in our case these constraints (27) are stated in the proof of Lemma 6). They may be some partial differential equations which increase the differential order of the energy. However, they also depend strongly on the geometry: for instance, we do not yet know whether it is possible to find a limit energy depending on a higher gradient of the displacement (third or higher order gradient material) by changing the distribution of the high rigidity inclusions.

For

  every Borel set D and u ∈ L 1 (D), we denote by |D| the Lebesgue measure of D and by -D udx the mean value of u on D: -D udx := |D| -1 D udx .

Fig. 1 .Fig. 2 .

 12 Fig. 1. The composite material, Ω

3 ( 1 y

 31 2 dy = 0 .

  . on B } is dense in H . Then, we can restrict our study to a function u ∈ H . As E 0 is continuous on H , the result can be generalized to H . Let us choose a sequence R ε such that r ε << R ε << ε, and let us divide M ε in two parts by introducing a transition layer C ε

  Thus the following estimation for ∂u ε /∂r can be derived Finally, it is easy to verify that |∂w p ε /∂x 3 | is bounded; then there exists M 4 such that , (cos θ, sin θ))| ≤ M 4 .(46)The estimations (44), (45), (46) imply that |∇u ε | is bounded on each layer C p ε , and thence on the set C ε . As |C ε | tends to 0, there follows lim sup Cε |∇u ε | 2 dx = 0 and lim sup E m

	|	∂u ε ∂r	(x p 1 + r cos θ, x p 2 + r sin θ, x 3 )|
					= |(1 -γ(r))	∂u ∂r	(x p 1 + r cos θ, x p 2 + r sin θ, x 3 )
					+	d γ dr	w p ε (x 3 , (cos θ, sin θ) -u(x p 1 + rcosθ, x p 2 + rsinθ, x 3 ) | ,
					≤ M + M 3 log(	r ε R ε	)	-1	.	(45)
	|	∂u ε ∂x 3	(x p 1 + r cos θ, x p 2 + r sin θ, x 3 )|
							= |(1 -γ(r))	∂u ∂x 3	(x p 1 + r cos θ, x p 2 + r sin θ, x 3 )
							+γ(r)	∂w p ε ∂x 3	(x 3

Acknowledgements. This paper was initiated by enlightening discussions with M. Bellieud and G. Bouchitté who studied non-local effects for scalar problems in the same geometry [5,4].