Franz Chouly
email: franz.chouly@univ-fcomte.fr

Patrick Hild
email: patrick.hild@math.univ-toulouse.fr

On convergence of the penalty method for unilateral contact problems

Keywords: unilateral contact, variational inequality, finite elements, penalty method, a priori error estimates AMS Subject Classification: 65N12, 65N30, 35J86, 74M15

We present a convergence analysis of the penalty method applied to unilateral contact problems in two and three space dimensions. We first consider, under various regularity assumptions on the exact solution to the unilateral contact problem, the convergence of the continuous penalty solution as the penalty parameter ε vanishes. Then, the analysis of the finite element discretized penalty method is carried out. Denoting by h the discretization parameter, we show that the error terms we consider give the same estimates as in the case of the constrained problem when the penalty parameter is such that ε = h.

Introduction

The penalty method is a classical and widespread method for the numerical treatment of constrained problems, in particular the unilateral contact problems arising in mechanics of deformable bodies which involve a nonlinear boundary condition written as an inequality (see, e.g., [START_REF] Kikuchi | Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF][START_REF] Laursen | Computational contact and impact mechanics[END_REF][START_REF] Wriggers | Computational Contact Mechanics[END_REF]). Nevertheless, and to the best of our knowledge, the convergence analysis of the method in the simplest case of linear elastostatics with or without finite element discretization has been object of few studies. We may nevertheless quote the earlier, and pioneering works of Kikuchi, Kim, Oden and Song [START_REF] Kikuchi | Penalty-finite-element approximation of a class of unilateral problems in linear elasticity[END_REF][START_REF] Oden | Finite element methods for constrained problems in elasticity[END_REF][START_REF] Oden | Interior penalty methods for finite element approximations of the Signorini problem in elastostatics[END_REF] (see also [START_REF] Kikuchi | Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF]) and the more recent study dealing with the boundary element method [START_REF] Chernov | A priori error estimates for hp penalty bem for contact problems in elasticity[END_REF]. In the context of elliptic partial differential equations, the penalty method is also classical for the treatment of Dirichlet boundary conditions, and has been thoroughly analyzed for instance in [START_REF] Babuška | The finite element method with penalty[END_REF][START_REF] Barrett | Finite element approximation of the Dirichlet problem using the boundary penalty method[END_REF]. However, due to the very different nature of Dirichlet boundary conditions and contact conditions, and of the resulting penalty methods, the aforementioned analysis can be hardly adapted to contact problems. In this paper, we present a convergence analysis of the penalty method for unilateral contact which use in particular some recent results from [START_REF] Hild | An improved a priori error analysis for finite element approximations of Signorini's problem[END_REF]. We analyze both the continuous and discrete problems. We limit the analysis to a conformal discretization with piecewise linear finite elements. In particular, we show that the same (sub-optimal, quasi-optimal or optimal) convergence rates as for the constrained problem (governed by a variational inequality) can be recovered with the choice ε = h, where ε is the penalty parameter, and h is the mesh size. A remarkable fact is that this choice is independent of the regularity of the continuous solution. Let us introduce first some useful notations. In what follows, bold letters like u, v, indicate vector or tensor valued quantities, while the capital ones (e.g., V, K, . . .) represent functional sets involving vector fields. As usual, we denote by (H s (•)) d , s ∈ R, d = 1, 2, 3 the Sobolev spaces in one, two or three space dimensions (see [START_REF] Adams | Sobolev spaces[END_REF]). The Sobolev norm of (H s (D)) d (dual norm if s < 0) is denoted by • s,D and we keep the same notation when d = 1, d = 2 or d = 3. The letter C stands for a generic constant, independent of the discretization parameters.

Setting

The contact problem

We consider an elastic body Ω in R d with d = 2 or d = 3. Small strain assumptions are made, as well as plane strain when d = 2. The boundary ∂Ω of Ω is polygonal or polyhedral and we suppose that ∂Ω consists in three nonoverlapping parts Γ D , Γ N and the contact boundary Γ C , with meas(Γ D ) > 0 and meas(Γ C ) > 0. The contact boundary is supposed to be a straight line segment when d = 2 or a polygon when d = 3 to simplify. The normal unit outward vector on ∂Ω is denoted n. In its initial stage, the body is in contact on Γ C with a rigid foundation and we suppose that the unknown final contact zone after deformation will be included into Γ C . The body is clamped on Γ D for the sake of simplicity. It is subjected to volume forces f ∈ (L 2 (Ω)) d and to surface loads g ∈ (L 2 (Γ N )) d . The unilateral contact problem in linear elasticity consists in finding the displacement field u : Ω → R d verifying the equations and conditions (1)-(3):

div σ(u) + f = 0 in Ω, σ(u) = A ε(u) in Ω, u = 0 on Γ D , σ(u)n = g on Γ N , (1)
where σ = (σ ij ), 1 ≤ i, j ≤ d, stands for the stress tensor field and div denotes the divergence operator of tensor valued functions. The notation ε(v) = (∇v+∇v T )/2 represents the linearized strain tensor field and A is the fourth order symmetric elasticity tensor having the usual uniform ellipticity and boundedness property. For any displacement field v and for any density of surface forces σ(v)n defined on ∂Ω we adopt the following decomposition in normal and tangential components:

v = v n n + v t and σ(v)n = σ n (v)n + σ t (v).
The nonlinear conditions describing unilateral contact on Γ C are:

u n ≤ 0, σ n (u) ≤ 0, σ n (u) u n = 0, (2)
and the frictionless condition is

σ t (u) = 0. ( 3
)
We introduce the Hilbert space:

V := v ∈ H 1 (Ω) d : v = 0 on Γ D .
The convex cone of admissible displacements which satisfy the noninterpenetration on the contact zone Γ C is:

K := {v ∈ V : v n = v • n ≤ 0 on Γ C } . Define a(u, v) := Ω σ(u) : ε(v) dΩ, L(v) := Ω f • v dΩ + Γ N g • v dΓ,
for any u and v in V. From the previous assumptions, we know that a(•, •) is bilinear, symmetric, V-elliptic and continuous on V × V. We know also that L(•) is a continuous linear form on V.

The weak formulation of Problem (1)-(3), as a variational inequality (see [START_REF] Fichera | Problemi elastostatici con vincoli unilaterali: Il problema di Signorini con ambigue condizioni al contorno[END_REF] and also [START_REF] Haslinger | Contact between elastic bodies. I. Continuous problems[END_REF][START_REF] Haslinger | Handbook of numerical analysis[END_REF][START_REF] Kikuchi | Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF]), is:

Find u ∈ K such that: a(u, v -u) ≥ L(v -u), ∀ v ∈ K. (4)
Stampacchia's Theorem ensures that Problem (4) admits a unique solution.

The penalty formulation of the unilateral contact problem

Let us first introduce the notation [•] + for the positive part of a scalar quantity a ∈ R:

[a] + = a if a > 0, 0 otherwise.
In the rest of this paper, we will make an extensive use of the following properties:

a ≤ [a] + , a[a] + = [a] 2 + , ∀a ∈ R. (5)
We introduce also:

W = v n |Γ C : v ∈ V ,
and its topological dual space W ′ , endowed with its usual dual norm. Since Γ C is a straight line segment or a polygon, we have

H 1/2 00 (Γ C ) ⊂ W ⊂ H 1/2 (Γ C ) which implies W ′ ⊂ (H 1/2 00 (Γ C )) ′ where H 1/2 00 (Γ C
) is the space of the restrictions on Γ C of functions in H 1/2 (∂Ω) vanishing outside Γ C and H 1/2 (Γ C ) is the space of the restrictions on Γ C of traces on ∂Ω of functions in H 1 (Ω). We refer to [START_REF] Lions | Non-homogeneous boundary value problems and applications[END_REF] and [START_REF] Adams | Sobolev spaces[END_REF] for a detailed presentation of trace operators and/or trace spaces. Let ε > 0 be a small parameter. The penalty method for unilateral contact problem (4) reads:

   Find u ε ∈ V such that: a(u ε , v) + 1 ε [u ε,n ] + , v n Γ C = L(v), ∀ v ∈ V, (6)
where •, • Γ C stands for the duality product between W ′ and W . Note that this formulation (which is an approximation of the exact contact conditions ( 2)) is obtained by setting the condition

σ n (u ε ) = -1 ε [u ε,n ] + instead of the conditions in (2) on the contact boundary Γ C . Remark 2.1. Note that for Problem (6), it holds : u ε,n , v n ∈ W ⊂ L 2 (Γ C ), and so [u ε,n ] + ∈ L 2 (Γ C ).
It results that the duality product •, • Γ C can also be understood simply as the inner product in L 2 (Γ C ). This remark still holds for the remaining part of the paper in which all the duality products •, • Γ C can be changed with inner products.

We recall that Problem ( 6) is well-posed using an argument proposed by H. Brezis for M type and pseudo-monotone operators [START_REF] Brezis | Équations et inéquations non linéaires dans les espaces vectoriels en dualité[END_REF] (see also [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF] and [START_REF] Kikuchi | Penalty-finite-element approximation of a class of unilateral problems in linear elasticity[END_REF]):

Theorem 2.2. For all ε > 0, Problem (6) admits a unique solution u ε .

Proof. Using the Riesz representation theorem, we define a (nonlinear) operator B : V → V with the following formula:

(Bv, w) 1,Ω := a(v, w) + 1 ε [v n ] + , w n Γ C , ∀v, w ∈ V,
where (., .) 1,Ω denotes the inner product in (H 1 (Ω)) d . Note that Problem ( 6) is well-posed if and only if B is a one-to-one operator. Let v, w ∈ V, it follows from the definition of B that:

(Bv -Bw, v -w) 1,Ω = a(v -w, v -w) + 1 ε [v n ] + -[w n ] + , v n -w n Γ C .
Due to the properties (5), we observe that, for all a, b ∈ R:

([a] + -[b] + )(a -b) = a[a] + + b[b] + -b[a] + -a[b] + ≥ [a] 2 + + [b] 2 + -2[a] + [b] + = ([a] + -[b] + ) 2 ≥ 0.
This property combined to the V-ellipticity of a(•, •) imply that there exists α > 0 such that:

(Bv -Bw, v -w) 1,Ω ≥ α v -w 2 1,Ω , ∀v, w ∈ V. (7)
Let us also show that the operator B is hemicontinuous, which means that for all v, w ∈ V, the real function

[0, 1] ∋ t → ϕ(t) := (B(v -tw), w) 1,Ω ∈ R is continuous. For s, t ∈ [0, 1],
we have:

|ϕ(t) -ϕ(s)| = |(B(v -tw) -B(v -sw), w) 1,Ω | ≤ |a(v -tw, w) -a(v -sw, w)| + 1 ε | [v n -tw n ] + , w n Γ C -[v n -sw n ] + , w n Γ C | = |a((s -t)w, w)| + 1 ε Γ C ([v n -tw n ] + -[v n -sw n ] + ) w n dΓ ≤ |t -s|a(w, w) + 1 ε Γ C |[v n -tw n ] + -[v n -sw n ] + | |w n | dΓ.
Using then the inequality

|[v n ] + -[w n ] + | ≤ |v n -w n |, we obtain: |ϕ(t) -ϕ(s)| ≤ |t -s|a(w, w) + 1 ε Γ C |t -s| |w n | |w n | dΓ = |t -s| a(w, w) + 1 ε w n 2 0,Γ C .
It follows that ϕ is Lipschitz, so continuous. The operator B is then hemicontinuous. Since (7) also holds, we can apply the Corollary 15 (p.126) of [START_REF] Brezis | Équations et inéquations non linéaires dans les espaces vectoriels en dualité[END_REF] to conclude that B is a one-to-one operator from V to V. That concludes the proof of the theorem.

Finite element setting and discrete penalty problem

Let V h ⊂ V be a family of finite dimensional vector spaces (see [START_REF] Ciarlet | Handbook of numerical analysis[END_REF]) indexed by h coming from a family T h of triangulations of the domain Ω (h = max T ∈T h h T where h T is the diameter of T ). The family of triangulations is supposed regular, i.e., there exists σ > 0 such that ∀T ∈ T h , h T /ρ T ≤ σ where ρ T denotes the radius of the inscribed ball in T . We choose standard continuous and piecewise affine functions, i.e.:

V h = v h ∈ (C(Ω)) d : v h |T ∈ (P 1 (T )) d , ∀T ∈ T h , v h = 0 on Γ D .
We introduce W h (Γ C ), the space of normal traces on Γ C for discrete functions in V h :

W h (Γ C ) := µ h ∈ C(Γ C ) : ∃v h ∈ V h , v h • n = µ h .
To simplify, we suppose that the end points (or the border for d = 3) of Γ C belong to Γ N . Moreover we assume that the mesh on Γ C induced by T h , on which are defined functions of W h (Γ C ), is quasi-uniform, which implies in particular that it is locally quasi-uniform, in the sense of Bramble & al [START_REF] Bramble | On the stability of the L 2 projection in H 1 (Ω)[END_REF].

The discrete version of the penalty method [START_REF] Bramble | On the stability of the L 2 projection in H 1 (Ω)[END_REF] for Problem (4) reads:

   Find u h ε ∈ V h such that: a(u h ε , v h ) + 1 ε [u h ε,n ] + , v h n Γ C = L(v h ), ∀ v h ∈ V h , (8)
where ε > 0 is still the small penalty parameter. Using exactly the same argument as in Theorem 2.2, we see that this problem admits one unique solution.

Convergence analysis of the penalty method

We present in this section the convergence analysis of both continuous and discrete penalty methods. We first state the two main theorems. Next, we give some necessary technical lemmas, followed by the proof of each theorem.

Main results

Convergence when ε → 0

The penalty formulation introduces a consistency error. The following theorem gives a bound for this error as a function of the penalty parameter ε. In particular, we recover the well-known result that the solution u ε of Problem ( 6) converges to the solution u of Problem (4) as ε → 0. We suppose that u is more than H 3 2 regular so that all the duality pairings become L 2 inner products.

Theorem 3.1. Suppose that u, the solution of Problem (4), belongs to (H 3 2 +ν (Ω)) d with ν ∈ (0, 1/2]. Let u ε be the solution of Problem [START_REF] Bramble | On the stability of the L 2 projection in H 1 (Ω)[END_REF]. We have the a-priori estimates:

u -u ε 1,Ω ≤ Cε 1 2 +ν u 3 2 +ν,Ω , (9)
σ n (u) + 1 ε [u ε,n ] + 0,Γ C ≤ Cε ν u 3 2 +ν,Ω , (10)
σ n (u) + 1 ε [u ε,n ] + -ν,Γ C ≤ Cε 2ν u 3 2 +ν,Ω , (11)
with C > 0 a constant, independent of ε and u.

3.1.2.

Convergence when ε → 0, h → 0 Next, the following theorem provides the convergence rates for the discrete penalty problem, as a function of both discretization parameters ε and h.

u -u h ε 1,Ω + ε 1 2 σ n (u) + 1 ε [u h ε,n ] + 0,Γ C ≤ C        h 1 2 + ν 2 +ν 2 + h ν ε 1 2 + h ν-1 2 ε u 3 2 +ν,Ω if 0 < ν < 1 2 ,
h| ln h|

1 2 + (hε) 1 2 + ε u 2,Ω if ν = 1 2 , (12)
σ n (u) + 1 ε [u h ε,n ] + -ν,Γ C ≤ C        h 1 2 + 3ν 2 +ν 2 ε -1 2 + h 3ν 2 +ν 2 + h 2ν-1 2 ε 1 2 + h 2ν-1 ε u 3 2 +ν,Ω if 0 < ν < 1 2 , h 3 2 | ln h| 1 2 ε -1 2 + h| ln h| 1 2 + (hε) 1 2 + ε u 2,Ω if ν = 1 2 , ( 13
)
with C > 0 a constant, independent of ε, h and u. In three space dimensions, the terms h

1 2 + ν 2 +ν 2 (resp. h| ln h| 1
2 ) in ( 12) and ( 13) have to be replaced with h 4 ).

1 2 + ν 2 (resp. h 3
One interesting fact is that the choice ε = h in the previous theorem, leads to the same error estimates as for the finite element approximation of the variational inequality (see [START_REF] Hild | An improved a priori error analysis for finite element approximations of Signorini's problem[END_REF]). With this choice, the error estimates are straightforward and are given next.

Corollary 3.3. Suppose that the solution u of Problem (4) belongs to (H 3 2 +ν (Ω)) d with ν ∈ (0, 1/2]. Suppose also that the penalty parameter is chosen as ε = h. The solution u h ε of the discrete penalty problem (8) satisfies the following error estimates in two space dimensions:

u -u h ε 1,Ω + h 1 2 σ n (u) + 1 ε [u h ε,n ] + 0,Γ C + h 1 2 -ν σ n (u) + 1 ε [u h ε,n ] + -ν,Γ C ≤ C      h 1 2 + ν 2 +ν 2 u 3 2 +ν,Ω if 0 < ν < 1 2 , h| ln h| 1 2 u 2,Ω if ν = 1 2 , ( 14
)
with C > 0 a constant, independent of h and u. In three space dimensions, the terms h

1 2 + ν 2 +ν 2 (resp. h| ln h| 1
2 ) in [START_REF] Haslinger | Handbook of numerical analysis[END_REF] have to be replaced with h 4 ).

1 2 + ν 2 (resp. h 3
Remark 3.4. Note that the quasi-optimality is only due to the estimation of the contact term, and is also present for standard finite element discretization [START_REF] Hild | An improved a priori error analysis for finite element approximations of Signorini's problem[END_REF]. This is not an intrinsic property of the penalty method. Note that optimality could be recovered in both 2D and 3D cases (and for any 0 < ν ≤ 1/2) if additional assumptions on the transition area between contact and non contact are made (see [START_REF] Hüeber | An optimal a priori error estimate for nonlinear multibody contact problems[END_REF]).

Remark 3.5. If u belongs to (H s (Ω)) 2 with s > 2 (Ω ⊂ R 2 ) then it is easy to show (by using the result in [START_REF] Renard | A quasi-optimal a priori error estimate for two-dimensional Signorini problem approximated by linear finite elements[END_REF]) that the error term in ( 14) is bounded by Ch u s,Ω .

Remark 3.6. To ensure convergence of the term u -

u h ε 1,Ω + ε 1 2 σ n (u) + 1 ε [u h ε,n ] + 0,Γ C ,
we need to choose at least ε < Ch 1 2 -ν (see [START_REF] Fichera | Problemi elastostatici con vincoli unilaterali: Il problema di Signorini con ambigue condizioni al contorno[END_REF]). As stated previously, to recover optimal convergence rate of the terms ( 12) and ( 13), the best choice is ε = h.

Preliminary technical lemmas

Let P h : L 2 (Γ C ) → W h (Γ C ) denote the L 2 (Γ C )-projection operator onto W h (Γ C ).
We recall in this lemma the stability and interpolation properties of P h : Lemma 3.7. Suppose that the mesh associated to W h (Γ C ) is locally quasi-uniform. For all s ∈ [0, 1] and all v ∈ H s (Γ C ), we have the stability estimate:

P h v s,Γ C ≤ C v s,Γ C . ( 15
)
The following interpolation estimate also holds:

v -P h v 0,Γ C ≤ Ch s v s,Γ C , (16)
for all v ∈ H s (Γ C ). The constant C > 0 is in both cases independent of v and h.

Proof. The stability estimate ( 15) is proven in [START_REF] Bramble | On the stability of the L 2 projection in H 1 (Ω)[END_REF]. The interpolation estimate comes from e.g. [START_REF] Bernardi | A new nonconforming approach to domain decomposition: the mortar element method[END_REF]. We need another lemma which concerns the existence of a discrete bounded lifting from Γ C to Ω:

Lemma 3.8. Suppose that the mesh on the contact boundary Γ C is quasi-uniform. There exists

R h : W h (Γ C ) → V h and C > 0, such that: R h (v h )| Γ C • n = v h , R h (v h ) 1,Ω ≤ C v h 1 2 ,Γ C , ( 17
)
for all v h ∈ W h (Γ C ).
Proof. The existence of such an operator is proven in [START_REF] Bjørstad | Iterative methods for the solution of elliptic problems on regions partitioned into substructures[END_REF] (see also [START_REF] Domínguez | Stability of discrete liftings[END_REF]). The proof of each theorem relies strongly on appropriate estimates for the approximate contact condition on Γ C in dual norm. We give these estimates in the following lemma: Lemma 3.9. Suppose that u, the solution of Problem (4), belongs to (H 3 2 +ν (Ω)) d with ν ∈ (0, 1/2]. Let u ε be the solution of Problem [START_REF] Bramble | On the stability of the L 2 projection in H 1 (Ω)[END_REF]. We have the bound:

σ n (u) + 1 ε [u ε,n ] + -ν,Γ C ≤ C ε ν σ n (u) + 1 ε [u ε,n ] + 0,Γ C + ε ν-1 2 u -u ε 1,Ω , (18)
with C > 0 a constant, independent of ε, u and u ε . Let u h ε be the solution of Problem [START_REF] Chernov | A priori error estimates for hp penalty bem for contact problems in elasticity[END_REF]. We have also the bound:

σ n (u) + 1 ε [u h ε,n ] + -ν,Γ C ≤ C h ν σ n (u) + 1 ε [u h ε,n ] + 0,Γ C + h ν-1 2 u -u h ε 1,Ω , (19)
with C > 0 a constant, independent of ε, h, u and u h ε .

Proof. We first prove the bound [START_REF] Laursen | Computational contact and impact mechanics[END_REF] for the discrete solution. The bound [START_REF] Kikuchi | Penalty-finite-element approximation of a class of unilateral problems in linear elasticity[END_REF] for the continuous solution will be derived in a quite similar manner. By definition of the norm • -ν,Γ C and using the projection operator P h :

σ n (u) + 1 ε [u h ε,n ] + -ν,Γ C = sup v∈H ν (Γ C ) σ n (u) + 1 ε [u h ε,n ] + , v Γ C v ν,Γ C ≤ sup v∈H ν (Γ C ) σ n (u) + 1 ε [u h ε,n ] + , v -P h v Γ C v ν,Γ C + sup v∈H ν (Γ C ) σ n (u) + 1 ε [u h ε,n ] + , P h v Γ C v ν,Γ C ≤ σ n (u) + 1 ε [u h ε,n ] + 0,Γ C sup v∈H ν (Γ C ) v -P h v 0,Γ C v ν,Γ C + C sup v∈H ν (Γ C ) σ n (u) + 1 ε [u h ε,n ] + , P h v Γ C P h v ν,Γ C ≤ C h ν σ n (u) + 1 ε [u h ε,n ] + 0,Γ C + sup v∈H ν (Γ C ) σ n (u) + 1 ε [u h ε,n ] + , P h v Γ C P h v ν,Γ C . (20)
In the fourth line, we used the Cauchy-Schwarz inequality and the stability property [START_REF] Hild | An improved a priori error analysis for finite element approximations of Signorini's problem[END_REF]. The last line is a direct consequence of the interpolation property [START_REF] Hüeber | An optimal a priori error estimate for nonlinear multibody contact problems[END_REF]. Now, with help of the relationship a(u

-u h ε , v h ) = σ n (u) + 1 ε [u h ε,n ] + , v h n Γ C for all v h ∈ V h ,
of the discrete lifting (17) and using the continuity of a(•, •), it results that:

sup v∈H ν (Γ C ) σ n (u) + 1 ε [u h ε,n ] + , P h v Γ C P h v ν,Γ C = sup v∈H ν (Γ C ) σ n (u) + 1 ε [u h ε,n ] + , R h (P h v)| Γ C • n Γ C P h v ν,Γ C = sup v∈H ν (Γ C ) a(u -u h ε , R h (P h v)) P h v ν,Γ C ≤ C u -u h ε 1,Ω sup v∈H ν (Γ C ) R h (P h v) 1,Ω P h v ν,Γ C ≤ C u -u h ε 1,Ω sup v∈H ν (Γ C ) P h v 1 2 ,Γ C P h v ν,Γ C .
Since the discrete trace space W h (Γ C ) has the quasi-uniform mesh property, we make use of the inverse inequality

P h v 1 2 ,Γ C ≤ Ch ν-1 2 P h v ν,Γ C ,
so that we finally get:

sup v∈H ν (Γ C ) σ n (u) + 1 ε [u h ε,n ] + , P h v Γ C P h v ν,Γ C ≤ Ch ν-1 2 u -u h ε 1,Ω .
This, together with (20) leads to:

σ n (u) + 1 ε [u h ε,n ] + -ν,Γ C ≤ C h ν σ n (u) + 1 ε [u h ε,n ] + 0,Γ C + h ν-1 2 u -u h ε 1,Ω ,
which is the desired bound [START_REF] Laursen | Computational contact and impact mechanics[END_REF].

For the continuous problem, we introduce V ε , a fictitious finite element space, defined identically as V h and with the choice of mesh size h = ε. We define also a fictitious discrete trace space W ε (Γ C ), in the same manner as W h (Γ C ), and with also the mesh size h = ε. We note simply

P ε : L 2 (Γ C ) → W ε (Γ C ) the L 2 (Γ C )-projection operator on W ε (Γ C ).
We again obtain an analogous estimate as in [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF] :

σ n (u) + 1 ε [u ε,n ] + -ν,Γ C ≤ C ε ν σ n (u) + 1 ε [u ε,n ] + 0,Γ C + sup v∈H ν (Γ C ) σ n (u) + 1 ε [u ε,n ] + , P ε v Γ C P ε v ν,Γ C .
With help of the relationship

a(u -u ε , v) = σ n (u) + 1 ε [u ε,n ] + , v n Γ C , v ∈ V, (21)
using the continuity of a(•, •) and a discrete lifting operator R ε , we come to the conclusion that:

σ n (u) + 1 ε [u ε,n ] + -ν,Γ C ≤ C ε ν σ n (u) + 1 ε [u ε,n ] + 0,Γ C + ε ν-1 2 u -u ε 1,Ω
which is the desired result [START_REF] Kikuchi | Penalty-finite-element approximation of a class of unilateral problems in linear elasticity[END_REF].

Remark 3.10. When ν = 1/2 the previous result can be improved, by using simply estimate [START_REF] Lions | Non-homogeneous boundary value problems and applications[END_REF] to obtain

σ n (u) + 1 ε [u ε,n ] + -1/2,Γ C ≤ C u -u ε 1,Ω .

Proof of Theorem 3.1

We first use the V-ellipticity of a(•, •) and the fact that u (resp. u ε ) is the solution of Problem (4) (resp. ( 6)) to obtain:

α u -u ε 2 1,Ω ≤ a(u -u ε , u -u ε ) = a(u, u -u ε ) -a(u ε , u -u ε ) = σ n (u) + 1 ε [u ε,n ] + , u n -u ε,n Γ C = σ n (u), u n Γ C + 1 ε [u ε,n ] + , u n Γ C -σ n (u), u ε,n Γ C - 1 ε [u ε,n ] + , u ε,n Γ C ,
where α > 0 denotes the ellipticity constant. Contact conditions (2) on Γ C yield:

σ n (u), u n Γ C = 0, 1 ε [u ε,n ] + , u n Γ C ≤ 0.
Once again conditions (2) associated to properties (5) provide:

-σ n (u), u ε,n Γ C ≤ -σ n (u), [u ε,n ] + Γ C , 1 ε [u ε,n ] + , u ε,n Γ C = 1 ε [u ε,n ] + , [u ε,n ] + Γ C .
We then get:

α u -u ε 2 1,Ω ≤ -σ n (u) + 1 ε [u ε,n ] + , [u ε,n ] + Γ C .
So we continue bounding:

α u -u ε 2 1,Ω ≤ -ε(σ n (u) + 1 ε [u ε,n ] + ), 1 ε [u ε,n ] + + σ n (u) -σ n (u) Γ C = -ε σ n (u) + 1 ε [u ε,n ] + 2 0,Γ C + ε σ n (u) + 1 ε [u ε,n ] + , σ n (u) Γ C ≤ -ε σ n (u) + 1 ε [u ε,n ] + 2 0,Γ C + ε δ σ n (u) + 1 ε [u ε,n ] + -ν,Γ C ε 1-δ σ n (u) ν,Γ C ≤ -ε σ n (u) + 1 ε [u ε,n ] + 2 0,Γ C + ε 2δ 2β σ n (u) + 1 ε [u ε,n ] + 2 -ν,Γ C + βε 2-2δ 2 σ n (u) 2 ν,Γ C , with δ ∈ [0, 1], β > 0. Note that since we supposed u ∈ (H 3 2 +ν (Ω)) 2 , we have σ n (u) ∈ H ν (Γ C
). We combine this result with estimation (18):

α u -u ε 2 1,Ω ≤ -ε 1 -C ε 2(δ+ν)-1 β σ n (u) + 1 ε [u ε,n ] + 2 0,Γ C + C ε 2(δ+ν)-1 β u -u ε 2 1,Ω + βε 2-2δ 2 σ n (u) 2 ν,Γ C ,
which can be transformed into:

α -C ε 2(δ+ν)-1 β u -u ε 2 1,Ω + ε 1 -C ε 2(δ+ν)-1 β σ n (u) + 1 ε [u ε,n ] + 2 0,Γ C ≤ βε 2-2δ 2 σ n (u) 2 ν,Γ C . Taking δ = 1/2 -ν, β = 2C max(1, α -1
) and using the estimate

σ n (u) ν,Γ C ≤ C u 3 2 +ν,Ω
proves the bounds ( 9) and ( 10) of the theorem. The bound (11) of the theorem is then a direct consequence of this last result and the estimation in norm of the negative exponent Sobolev space (18).

Proof of Theorem 3.2

We denote by I h the Lagrange interpolation operator mapping onto V h . We first use the V-ellipticity and the continuity of a(•, •), as well as Young's inequality, to obtain:

α u -u h ε 2 1,Ω ≤ a(u -u h ε , u -u h ε ) = a(u -u h ε , (u -I h u) + (I h u -u h ε )) ≤ C u -u h ε 1,Ω u -I h u 1,Ω + a(u -u h ε , I h u -u h ε ) ≤ α 2 u -u h ε 2 1,Ω + C 2 2α u -I h u 2 1,Ω + a(u, I h u -u h ε ) -a(u h ε , I h u -u h ε ),
with α > 0 the ellipticity constant. Since u is solution of (4) and u h ε is solution of (8), we can transform the term a(u,

I h u -u h ε ) -a(u h ε , I h u -u h ε ). So we obtain: α 2 u -u h ε 2 1,Ω ≤ C 2 2α u -I h u 2 1,Ω + σ n (u) + 1 ε [u h ε,n ] + , (I h u) n -u h ε,n Γ C .
Because of conditions (3) and since the 1D-Lagrange interpolation with piecewise-linear polynomials preserves the positivity, we have that (I h u) n ≤ 0 on Γ C . This implies:

1 ε [u h ε,n ] + , (I h u) n Γ C ≤ 0.
Once again condition (3) associated to properties (5) yield:

-σ n (u), u h ε,n Γ C ≤ -σ n (u), [u h ε,n ] + Γ C , 1 ε [u h ε,n ] + , u h ε,n Γ C = 1 ε [u h ε,n ] + , [u h ε,n ] + Γ C .
This results into:

α 2 u -u h ε 2 1,Ω ≤ C 2 2α u -I h u 2 1,Ω + σ n (u), (I h u) n Γ C -σ n (u) + 1 ε [u h ε,n ] + , [u h ε,n ] + Γ C . ( 22
)
We bound the last term σ n (u) + 1 ε [u h ε,n ] + , [u h ε,n ] + Γ C as follows:

-

σ n (u) + 1 ε [u h ε,n ] + , [u h ε,n ] + Γ C = -ε(σ n (u) + 1 ε [u h ε,n ] + ), 1 ε [u h ε,n ] + + σ n (u) -σ n (u) Γ C = -ε σ n (u) + 1 ε [u h ε,n ] + 2 0,Γ C + ε σ n (u) + 1 ε [u h ε,n ] + , σ n (u) Γ C ≤ -ε σ n (u) + 1 ε [u h ε,n ] + 2 0,Γ C + ε σ n (u) + 1 ε [u h ε,n ] + -ν,Γ C σ n (u) ν,Γ C .
The last inequality has been made possible due to the assumption on the regularity of u, which implies that σ n (u) ∈ H ν (Γ C ). Treating the last term with Young's inequality, and with β > 0, we obtain finally:

-σ n (u) + 1 ε [u h ε,n ] + , [u h ε,n ] + Γ C ≤ -ε σ n (u) + 1 ε [u h ε,n ] + 2 0,Γ C + ε 2 2β σ n (u) + 1 ε [u h ε,n ] + 2 -ν,Γ C + β 2 σ n (u) 2 ν,Γ C . (23)
We now combine this last inequality [START_REF] Oden | Interior penalty methods for finite element approximations of the Signorini problem in elastostatics[END_REF] with [START_REF] Oden | Finite element methods for constrained problems in elasticity[END_REF] and then insert the estimation [START_REF] Laursen | Computational contact and impact mechanics[END_REF] into the resulting inequality:

α 2 u -u h ε 2 1,Ω ≤ C 2 2α u -I h u 2 1,Ω + σ n (u), (I h u) n Γ C -ε σ n (u) + 1 ε [u h ε,n ] + 2 0,Γ C + ε 2 2β σ n (u) + 1 ε [u h ε,n ] + 2 -ν,Γ C + β 2 σ n (u) 2 ν,Γ C ≤ C 2 2α u -I h u 2 1,Ω + σ n (u), (I h u) n Γ C -ε 1 -C εh 2ν 2β σ n (u) + 1 ε [u h ε,n ] + 2 0,Γ C + Ch 2ν-1 ε 2 2β u -u ε 2 1,Ω + β 2 σ n (u) 2 ν,Γ C .
We rearrange the terms:

α 2 -Ch 2ν-1 ε 2 2β u -u h ε 2 1,Ω + ε 1 -C εh 2ν 2β σ n (u) + 1 ε [u h ε,n ] + 2 0,Γ C ≤ C 2 2α u -I h u 2 1,Ω + σ n (u), (I h u) n Γ C + β 2 σ n (u) 2 ν,Γ C . (24)
We choose β as follows:

β = C h 2ν-1 ε 2 α -1 + h 2ν ε ,
with C > 0 a constant sufficiently large so that the two left terms in [START_REF] Renard | A quasi-optimal a priori error estimate for two-dimensional Signorini problem approximated by linear finite elements[END_REF] are positive irrespectively of the values of ε and h. It results that:

β σ n (u) 2 ν,Γ C = C h 2ν-1 ε 2 + h 2ν ε σ n (u) 2 ν,Γ C . ( 25
)
The estimation of the Lagrange interpolation error in L 2 and H 1 norms on a domain D is classical (see e.g., [START_REF] Ern | Theory and practice of finite elements[END_REF]):

Theorem 3 . 2 . 3 2

 323 Suppose that the solution u of Problem (4) belongs to (H +ν (Ω)) d with ν ∈ (0, 1/2]. The solution u h ε of the discrete penalty problem (8) satisfies the following error estimates in two space dimensions:

for s ∈ [START_REF] Adams | Sobolev spaces[END_REF][START_REF] Babuška | The finite element method with penalty[END_REF]. The contact term σ n (u), (I h u) n Γ C can be estimated in two space dimensions using results from [START_REF] Hild | An improved a priori error analysis for finite element approximations of Signorini's problem[END_REF]:

In three space dimensions the bound is obtained in a straightforward way using (26) for any 0 < ν ≤ 1/2 :

In two space dimensions, we combine finally the estimations ( 24)-( 27) to prove that:

which is the first required estimate. This, together with the estimate [START_REF] Laursen | Computational contact and impact mechanics[END_REF], yields additionally the bound in two space dimensions:

(h

Using h 2ν ≤ h 3ν 2 +ν 2 and h ≤ h| ln h| 1 2 ends the proof of ( 13). The bounds in three space dimensions are obtained as before using estimate (28) instead of (27).


François Fouquet

Erwan Daubert

Noël Plouzeau

Olivier Barais

Johann Bourcier

Jean-Marc Jézéquel

Dissemination of reconfiguration policies on mesh networks

Component-based platforms are widely used to develop and deploy distributed pervasive system that exhibit a high degree of dynamicity, concurrency, distribution, heterogeneity, and volatility. This paper deals with the problem of ensuring safe yet efficient dynamic adaptation in a distributed and volatile environment. Most current platforms provide capabilities for dynamic local adaptation to adapt these systems to their evolving execution context, but are still limited in their ability to handle distributed adaptations. Thus, a remaining challenge is to safely propagate reconfiguration policies of component-based systems to ensure consistency of the architecture configuration models over a dynamic and distributed system. In this paper we implement a specific algorithm relying on the models at runtime paradigm to manage platform independent models of the current system architecture and its deployed configuration, and to propagate reconfiguration policies. We evaluate a combination of gossip-based algorithms and vector clock techniques that are able to propagate these policies safely in order to preserve consistency of architecture configuration models among all computation nodes of the system. This evaluation is done with a test-bed system running on a large size grid network.

Introduction

Nowadays, the increasing use of Internet of Things devices for computer supported cooperative work leads to large systems. As these devices use multiple mobile networks, these systems must deal with concurrency, distribution, and volatility issues. This volatility requires dynamic auto-adaptation of the system architecture, in order to provide domain specific services continuously. Tactical information and decision support systems for on field emergency management are perfect examples of such highly dynamic systems. Indeed, these multi-user interactive systems built on mobile devices need frequent changes of architecture to deal with rapid system evolution (e.g. scale up or scale down of team, download of new software modules by the device user) or to cope with network disconnections. For such systems, the traditional design process "design, code, compile, test, deploy, use, iterate" does not work.

Dynamic adaptation, pursuing IBM's vision of autonomic computing, is a very active area since the late 1990's -early 2000's [START_REF] Kephart | The Vision of Autonomic Computing[END_REF]. Modern componentbased systems [START_REF] Raj | Implementing service-oriented architecture (soa) with the java ee 5 sdk. Sun Microsystems[END_REF][START_REF] Bruneton | The fractal component model and its support in java: Experiences with auto-adaptive and reconfigurable systems[END_REF] provide a reflection and intercession layer to dynamically reconfigure a running system. But the reconfiguration process remains complex, unreliable and often irreversible in a volatile and distributed context. The use of model-driven techniques for managing such run-time behavior (named models@runtime [START_REF] Blair | Models@runtime[END_REF]) helps to handle software reconfiguration. Models@runtime basically pushes the idea of reflection [START_REF] Morin | Models@ run.time to support dynamic adaptation[END_REF] one step further by considering the reflection layer as a real model that can be uncoupled from the running architecture (e.g. for reasoning, validation, and simulation purposes) and later automatically resynchronized with its running instance to trigger reconfigurations. Kevoree is our open-source dynamic component model 1 , which relies on models at runtime to properly support the dynamic reconfiguration of distributed systems. The model used at runtime reflects the global running architecture and the distributed topology. In Kevoree, when a distributed node receives a model update that reflects the target running architecture, the node extracts the reconfigurations that affect it and transform them into a set of platform reconfiguration primitives. Finally, it executes them and propagates the reflection model to other nodes as a new consistent architecture model.

In a highly distributed and volatile environment, one of the challenges is the propagation of reconfiguration policies. Handling concurrent updates of shared data is a second challenge to be solved, as two nodes can trigger concurrent reconfigurations. Consistent dissemination of models at runtime in distributed systems requires a synchronization layer that solves these two challenges: information dissemination and concurrent update. Research in the field of peerto-peer communication has produced many algorithms to deal with information dissemination in a volatile context [START_REF] Eugster | From epidemics to distributed computing[END_REF]. Many paradigms are available to deal with this concurrent data exchange problems (e.g. vector clocks [START_REF] Fidge | Timestamps in message-passing systems that preserve the partial ordering[END_REF]). In this paper, we adapt a combination of gossip-based algorithms and vector clocks techniques to safely propagate reconfiguration policies by preserving architecture models consistency between all computation nodes of a distributed system. We have implemented a specific algorithm, which propagates configuration changes in a consistent manner in spite of frequent node link failures, relying on its payload of configuration data to improve its efficiency. We provide qualitative and quantitative evaluations of this algorithm, to help answering the following questions: (i) What is the influence of communication strategy on the propagation delay of models? (ii) Does a high rate of node link failure prevent the propagation of models and what is the impact of link failures on propagation delays? (iii) Does the algorithm detect concurrent updates of models and does it handle reconciliation correctly?

The remainder of this paper is organized as follows. Section 2 presents the background of this work. Section 3 details the combination of a gossip-based algorithm and the vector clock techniques used to preserve architecture models consistency between all computation nodes of the system. Section 4 details our experiments to evaluate this combination. Section 5 discusses articles, ideas and experimental results related to our work. Finally, Section 6 concludes this paper and presents ongoing work.

Background

Kevoree is an open-source dynamic component model 1 , which relies on models at runtime [START_REF] Blair | Models@runtime[END_REF] to properly support the dynamic adaptation of distributed sys- tems. Figure 1 presents a general overview of models@runtime. When changes appear as a new model (a target model) to apply on the system, it is checked and validated to ensure a well-formed system configuration. Then it will be compared with the current model that represents the running system. This comparison generates an adaptation model that contains the set of abstract primitives to go from the current model to the target one. Finally, the adaptation engine executes configuration actions to apply these abstract primitives. If an action fails, the adaptation engine rollbacks the configuration to ensure system consistency. Kevoree has been influenced by previous work that we carried out in the DiVA project [START_REF] Morin | Models@ run.time to support dynamic adaptation[END_REF]. With Kevoree we push our vision of models@runtime [START_REF] Morin | Models@ run.time to support dynamic adaptation[END_REF] farther.

In particular, Kevoree supports distributed models@runtime properly. To this aim we introduce the Node concept in the model to represent the infrastructure topology. Kevoree includes a Channel concept to allow for multiple communication semantics between remote Components deployed on heterogeneous nodes. All Kevoree concepts (Component, Channel, Node) obey the Type Object pattern [START_REF] Johnson | The Type Object Pattern[END_REF] to separate deployment artifacts from running artifacts. Kevoree supports multiple kinds of execution node technology (e.g. Java, Android, Mini-Cloud, FreeBSD, Arduino, . . . ). Kevoree also introduces a dedicated concept named Group, to encapsulate platform synchronization algorithms. Group allows to define communication channels between nodes to propagate reconfiguration policies (i.e. new target model). This Group concept also encapsulates a dedicated protocol to ensure specific synchronization policies (e.g. Paxos derived algorithms for total order synchronization, Gossip derived algorithms for partial order and opportunistic synchronization). Groups can be bound to several nodes (named members), allowing them to explicitly define different synchronization strategies for the overall distributed system. This architecture organization is illustrated in Figure 2. In addition, a Group also defines a scope of synchronization, i.e. it defines which elements of the global model must be synchronized for the group's members. This avoids to globally share model@runtime models.

P2P algorithm and mesh network Schollmeier [START_REF] Schollmeier | A definition of peer-to-peer networking for the classification of peer-to-peer architectures and applications[END_REF] defines a peer-to-peer network as "a distributed network architecture, where participants of the network share a part of their resources, which are accessible by the other peers directly, without passing intermediary entities". He also provides the following distinction: hybrid peer-to-peer networks use a central entity, while pure peerto-peer networks have no such entity. According to Wikipedia, a mesh network is "a type of network where each node must not only capture and disseminate its own data, but also serve as a relay for other nodes, that is, it must collaborate to propagate the data in the network". In these network topologies, gossip-like algorithms are good solutions to disseminate data.

Concurency data management for distributed message passing applications Distributed systems consist of a set of processes that cooperate to achieve a common goal. Processes communicate with data exchanges over the network, with no shared global memory. This leads to well known and difficult problems of causality and ordering of data exchanges. Solutions are known to cope with this problem: Lamport [START_REF] Lamport | Time, clocks, and the ordering of events in a distributed system[END_REF] defines an event order using logical clocks by adding a logical time on each message sent. Another solution was coinvented by Fidge [START_REF] Fidge | Timestamps in message-passing systems that preserve the partial ordering[END_REF] and Mattern [START_REF] Mattern | Virtual time and global states of distributed systems[END_REF], using a vector of logical clocks. In many cases the vector clock technique is the most appropriate solution to manage a partial order and concurrency between events [START_REF] Baldoni | Fundamentals of distributed computing[END_REF], e.g. in distributed hash table systems such as Voldemort2 ).

Synthesis In our vision of distributed environments, system management is decentralized, allowing each peer to build, maintain and alter the overall architecture and platform models at runtime. Because of nodes volatility, ensuring consistency during reconfiguration is a critical task. We use Kevoree and the notion of Group to encapsulate platform synchronization algorithms with gossip and vector clock techniques.

An algorithm to disseminate reconfiguration policies

Each node holds a copy of the model that describes the overall system configuration. This system model contains a description of the nodes that currently compose the system, of components that are installed on each node and of network links between nodes. It also contains all information about groups. A group is the unit of model consistency for the models at runtime technique. Each node involved in model consistency includes several named group instances, which participates in the distributed model management for the local node. Part 1 of the algorithm provides the data definition for one node.

In addition to the information given by the model, each group instance maintains specific information (see algorithm's Part 2): a group id, a local copy of the model and the local node id. It also stores its current vector clock, a score for each of its neighbors and a boolean attribute to record whether the model has changed since the last time another node requested the local node's vector clock. The score of the neighbors is used to select the more interesting one when the local node looks for new reconfigurations. When a change appears on the model stored in a node, the corresponding group instance is notified. The group instance then sends notification to all its neighbors. These neighbors in turn may send a message to the current node, to ask for model update information. As the underlying communication network is volatile and unreliable, some notifications can be lost and not received by some members of a group. To deal with these losses, each member of a group asks periodically a chosen group member for changes. Since a model is a rather large data, group instances ask for the vector clock of the remote instance first, in order to decide if a model transfer is needed. More precisely, after comparing the vector clock received with its own vector clock, a group instance will request a model if both vector clocks are concurrent or if the vector clock received is more recent than its local one. Here concurrency means that each local and remote model have different changes which dot not appear on the other. A vector clock is more recent than another if some changes appear on it but not on the other. Upon reception of a model, the group instance compares the model's vector clock and the local clock again. If the local vector clock is older, the local node updates its local clock and also updates the local copy of the model using the model just received. If the vector clocks are concurrent then the group must resolve this concurrency at the model level to compute the correct model and then update the vector clock accordingly.

Functions SelectPeer (see Algorithm Part 4) In addition to this mechanism, each node pulls periodically one of its neighbors, in order to cope for lost notifications. The selection of the neighbor to pull is controlled by a score mechanism: a score is assigned to each peer by the group instance and the selection of the peer is done according to the smaller score. The score of the node grows when it is selected or when the network link to access this node seems to be down. down link detection relies on a synchronization layer. This layer uses model information to check for all available peers periodically and then to notify the group instance of unreachable nodes. A peer score takes into account the duration of unavailability of the peer. When the peer becomes available, this number is reset to 0: restored availability clears the failure record. Indeed, as the system uses a sporadic and volatile network, peers often appear and disappear and most of the time disappearance events are not causally connected.

Functions about vector clocks (see Algorithm Part 5) Our algorithm relies on vector clocks to detect changes in remote configuration models. When a local update of the model appears, a boolean called changed is set to true to ensure that upon a vector clock request from another node the group instance will increment by one its version id in its local vector clock before sending it to the requesting peer. In case of concurrent updates of models we rely on the use of the reflexivity provided by the model at runtime to solve the conflict. Priority is given to information about the nodes already reached and affected by the update. Any node detecting a conflict will merge these models and their associated vector clocks to store it as its current state. A reasoning upper layer will then compute an update from this merged model by reading the model and correcting it. Description of this reasoning layer is beyond the scope of this paper and vector clocks merge and comparison is already defined on previous works on vector clocks [START_REF] Fidge | Timestamps in message-passing systems that preserve the partial ordering[END_REF] and [START_REF] Mattern | Virtual time and global states of distributed systems[END_REF].

Algorithm

Evaluation

We have performed qualitative and quantitative evaluations of our algorithm, aiming at measuring the following indicators: (1) model propagation delay;

(2) resilience to node link failure; (3) ability to detect concurrent models and to handle reconciliation. For each indicator we have set up an experimental protocol, using the firefighter tactical information case study metrics to simulate the system behaviour on a grid in different configurations. Although target platforms will be pervasive embedded systems, we have chosen a large scale grid as an evaluation testbed. The use of a grid allows us to stress the algorithm by setting up a large number of nodes but it also brings us more control over the parameters of the experiment e.g. network failure simulations. In this way experiments are reproducible, and reproducibility is essential to our experimental protocol. On-field validation is an ongoing work.

Common experimental protocol

Validation experiments share a common experimental protocol. Each experiment uses a set of logical Kevoree node deployed on physical nodes within a computer grid. Each Kevoree logical node is instantiated in a separate Java Virtual Machine and use the reference Kevoree implementation for JavaSE. The experimental grid is an heterogeneous grid that contains nodes of mixed computational power and type. Each node is connected to a local area network at 100 MB/s. Topology model All our experiments take a bootstrap model as input, which describes the current abstract architecture (i.e. in its platform independent form). This abstract model contains information on node representations, node logical links and node communication group instances and relationships. This node set and these relationships describe a topology of the system, which is used by our synchronization algorithm. In order to improve the simulation of a firefighter tactical information case study, we use a random generator to create topology models that are organized in a cluster of clusters. In this way it is easier to simulate non-direct communication (i.e. node A cannot communicate directly with node B but must pass through node C).

Global time axis traces

In order to track the propagation of new configurations in this distributed system, we decorate the algorithm with a logger. This logger sends a trace for each internal state change (i.e. new configuration or configuration reconciliation). These traces describe the current state of the group, namely the new vector clock, the identification of the peer originator of change and the network metrics used. In order to exploit temporal data on these traces without ensuring a global grid time synchronization we use a logger with a global time axis based on Java Greg Logger 3 . More precisely, this type of logger is based on a client server architecture. The server clock manages the global time reference. All clients periodically synchronize with the server, allowing it to store client latencies by taking into account clocks shift and network time transfer observed. Traces are emitted asynchronously by the client to the server, which then makes time reconciliation by adding the last value of latency observed for this client. All traces emitted by the server are therefore all time stamped accurately with the clock of the server. Finally, traces are chained by an algorithm to meet the following heuristic: a trace follows another one if it is the first occurrence that contains in its vector clock the originator node with its precise version number. Thus the final result for each experiment is a linked trace list on which we can precisely compute temporal results.

Communication modes We reuse mainly two classical exchange patterns to build our algorithm.

Pooling period term is associated with the time elapsed between two active synchronizations, and is initiated by a group member to another. In this synchronization step a vector clock and/or a model is sent back to the initiator. The Push/Pull technique is an association of the pooling active synchronization and an event-driven notification mechanism. This operation adds to the pooling mode a sending step to every reachable group member.

Experimental studies

Propagation delay versus network usage This first experiment aims at performing precise measurements of the capacity to disseminate model configurations. These measures will take care of the propagation delay and the network usage properties.

Experimental protocol As described in the common protocol subsection, measurements are performed on a computer grid. The probes injected in the Java implementation collect propagation delay and network occupation. After a bootstrap step on a topology model, a node chosen at random reconfigures its local model with a simple modification. In practice this reconfiguration step computes a new model, moving a component instance from one node to another chosen randomly. This new model is stored in the node, and the reconfiguration awaits propagation by the algorithm. This new configuration is tagged with the identification of reconfiguration originator. Figure 3 shows the topology model used for multi hop communication in the 66 nodes of this configuration. In this experiment, the network topology is static. No node joins or leaves the system. The experiment is driven by the following parameters:(1) delay before starting an active check of the peers update (model synchronization);(2) activation of sending of change notification messages. To evaluate the impact of the second parameter, the experiment is run twice. In the first run, notifications are not used and the active synchronization delay is set to 1000 ms. In the second run, notifications are used and active synchronization delay is 15 s. In both cases, a reconfiguration is triggered every 20 seconds and each reconfiguration run takes 240 seconds, resulting in 12 reconfiguration steps.

Analysis The observed per hop propagations delays are presented as a percentile distribution (see Graph 4). The values displayed are the raw values of absolute time logged divided by the minimum number of hops between the target and of the reconfiguration (the minimum being computed using a Bellman-Ford algorithm [START_REF] Cheng | A loopfree extended bellman-ford routing protocol without bouncing effect[END_REF]). The traffic volume from protocol messages is shown in Figure 5 in KB per node per reconfiguration; the volume does not include payload. Absolute values of network consumption depends highly of implementation. Results presented here are from the Java version and can be vastly improved when targeting embedded devices like microcontrollers.

The use of notification reduces the propagation delay significantly: the average value decreases from 1510 ms/hop to 215 ms/hop. In addition, percentile distribution shows that the standard deviations of propagations are lower with in the version with notification. Thus this version of the algorithm has better scalability for large graph diameters.

However, in comparing the push pull and the push algorithm, the use of notification on network usage is not as significant. Analysis shows that these results are affected by cycles in the topology. When using notification of change, nodes in cycles will create parallel branches configuration diffusion. This in turn will increase the number of conflict resolution to be done, and these resolutions increase network consumption unnecessarily, by exchanging the same model version. When notifications are not used, pooling delays are large enough to avoid this concurrent configuration "flood". As the payload is a model with topology information, the notification algorithm could use this information to prevent flood. This solution will be studied in future work.

Failures impact on propagation delay A mobile mesh network such as the one used in a firefighter tactical information system is characterized by a large number of nodes that are often unreachable. We designed our algorithm to cope with these network problems. The second experiment described below tests the ability of the algorithm to disseminate new models in a mesh network with different failure rates.

Experiment protocol

The experiment protocol is similar to the first experiment's one. The topology model is enhanced to provide a mesh network with many different routes between nodes (see Figure 6). At each run a modified model is pushed on a random node. The reconfiguration is similar to the previous experiment. During each run, additional failures are simulated on links between two nodes, according to a Poisson distribution. The failure rate is increased at each run, thus the number of reachable nodes decreases. To perform this failure simulation we inject probes, which also monitor synchronization events. At each run, the list of theoretically reachable nodes is computed and the initiator node waits for synchronization events from theses nodes. When all events have been received we compute the average propagation delay. In short, this experiment aims at checking that every theoretically reachable node receives the new configuration.

Analysis Figure 7 shows results of experiment #2. The histogram shows the rate of network failure for each run. The red curve displays the average propagation delay to reachable nodes (in milliseconds). Above a network failure of 85% the node originator of the reconfiguration is isolated from the network and therefore we stop the execution. With a failure rate under 85% every node receives the new configuration and we can compute the propagation delay.

Concurrency reconfiguration reconciliation Our third experiment addresses the problem of reconciliation and conflict detection between concurrent model updates. This problem occurs often in the firefighter tactical information case study architecture because of the sporadic communication capabilities of our network of nodes. As a node can stay isolated for some time, reconfiguration data no longer reaches it. Furthermore, local reconfigurations can also occur in its subnetwork. Connection restoration may produce conflicting concurrent model updates. We rely on vector clocks to detect these conflicts and on the conflicting model updates themselves. Experiment #3 aims at checking the behaviour of our algorithm in this conflicting updates situation.

Experiment protocol

The experiment protocol is based on experiment #2. We use a similar grid architecture but with only 12 nodes. An initial reconfiguration (c1) is launched on the p00 node just after the bootstrap phase. All network links are up. Then a fault is simulated on the link between nodes p00 and o00. Nodes o00, o01, o02 are then isolated. A new model is then pushed on node p00 (c2) and a different one on node o00 (c3). A delay of 1000 ms separates each reconfiguration and the algorithm is configured with a notification and a pooling period of 2000 ms.

Analysis Figure 8 shows results of experiment #3, which are derived from our branching algorithm traces. Three reconfigurations are represented as a succession of segments that show the propagation of updates. The first reconfiguration on the healthy network is represented in black (at time 0). Reconfiguration pushed on o00 (at time 2500) is represented in blue and the second reconfiguration pushed on p00 (time 2000) in red. The first reconfiguration propagates seamlessly to all nodes. At time 1500 a network failure is simulated. The second model given to p00 is propagated to all nodes except nodes reachable through o00 only. Similarly, the second model pushed on node o00 is not propagated to nodes after p00. At time 8000 we cancel the network failure simulated at time 1500. After a synchronization delay (380ms) we observe the branching of the two concurrent models as well as propagation of the merged version (purple line).

Discussion and related work

Our approach is dedicated to model at runtime synchronization, and combines commonly used paradigms in distributed computing like vector clocks (e.g. used in distributed hash table frameworks) and gossiping (e.g. used in social network graph dissemination). This section discusses our experimental results and compare them to other related work.

Vector clock size. Our first experiment measures the size of data exchanged during the reconfiguration step, as well as the time required to perform this reconfiguration. Figure 5 shows that the model@runtime synchronization overhead is significant, and this is mostly due to vector clock size. Many studies aim at reducing the data size of vector clocks, especially when synchronizing an unbounded number of peers. Sergio and al [START_REF] Almeida | Interval tree clocks: A logical clock for dynamic systems[END_REF] proposed the Interval Tree Clocks to optimize the mapping between the node identifier and its version. Our algorithm takes advantage of the model payload to allocate dynamic identifiers to nodes. Data such as node names or network identifications are stored in the payload itself and with this information we can already improve vector clocks. However, we plan to implement the interval tree clocks' fork and join model in the future. The size of exchanged data depends on the number of nodes and therefore modularization techniques are needed to maintain scalability and manage large mesh networks. Our approach addresses this need by exploiting the group structure of Kevoree. Each group instance synchronizes with a subset of nodes only, to keep the size of the vector clock under control.

Distributed reconfiguration capability. Concurrency management is a key problem in distributed systems. Many peer to peer systems solve it by having a single point of update for a given piece of data, limiting concurrent access to a one writer/many readers situation for that data. Realistic distributed configuration management is a many writers/many readers situation, because reconfigurations often involve more than one node. The simplest solution to this problem would use a single point for new configuration computation and dissemination start. As it avoids concurrency, such a system has a central point of failure incompatible with our use case. More advanced approaches such as the one proposed in [START_REF] Sykes | Flashmob: distributed adaptive self-assembly[END_REF] use distributed coordination techniques such as consensus to build the new configuration. They proposed an approach that allows the distributed nodes to collaborate to build the new configuration. Each node is responsible for building its local configuration. Configuration propagation is then done using a gossip algorithm without the need of concurrency management, since new configurations can be disseminated from a single originator node only. This approach based on a single source is unusable in our use case, because the sporadic nature of the nodes prevent their participation in a global consensus. On the contrary, our technique presented in this paper lets the distributed configuration evolve freely, even for nodes are isolated in unreachable groups. Every node can then compute a new global model that can be issued concurrently. Some approaches in distributed hash table implementations also rely on fully distributed data dissemination, e.g. Voldemort, where table modifications can occur in several nodes. This allows for service operation in degraded mode in the case of node disconnections. However, concurrency management must be managed separately. GossipKit [START_REF] Lin | Facilitating gossip programming with the gossipkit framework[END_REF] proposes a generic framework to evaluate and simulate gossip-derived algorithms. The project contains a minimal extensible event-based runtime with a set of components to define dedicated gossip protocols. We plan to integrate the GossipKit API in order to evaluate our algorithm on a GossipKit simulator.

Inverted communication and propagation delay. In our approach we reverse the traditional communication strategy of a gossip algorithm (push approach). New configurations are not directly pushed to the neighbours but they are stored instead, waiting for an active synchronisation by the neighbour (pull approach). This strategy lessens the impact of down network links on propagation delay, as shown by our experiment results on Figure 7. In addition, this enables message replay because a configuration is stored until neighbor connectivity is reestablished. These two properties are particularly useful for unreliable mesh networks. However, pull approaches have higher propagation time, but when combined with an observer pattern (a lazy push/pull approach) our results show that the gains are significant while keeping the interesting properties of pull. This experimental result is consistent with Leitao et al [START_REF] Leitão | Gossip-based broadcast[END_REF], which details several communication strategy for gossip algorithms.

Conclusion

In this paper we proposed a peer to peer and distributed dissemination algorithm to manage dynamic architectures based on the models at runtime paradigm. This algorithm is part of a larger framework that manages the continuous adaptation of pervasive systems. Using experimental results we have shown how our approach enhances reliability and guarantee of information delivery, by mixing and specializing different distributed algorithms. Our propagation algorithm relies on its payload (a model of the system) to overcome limits of vector clocks and to handle peer to peer concurrency conflicts. Thanks to the protocol layer based on vector clocks, a system architecture model propagated by the algorithm is always consistent, even on complex mesh network topologies. When concurrent updates are detected, the model at runtime layer is able to reconcile these updates to provide a valid architecture. By allowing each node to compute a new configuration, our approach supports dynamic adaptation on peer to peer networks without any central point of failure. This experimental demonstration of resilience on sporadic networks allows integration of our approach into adaptive architectures such as a firefighters tactical information system. In this direction, we are currently designing a dynamically scalable tactical information system in collaboration with a department of firefighters of Brittany; this system is a multi-user, real time decision system for incident management. 4 .

1Fig. 1 .Fig. 2 .

 12 Fig. 1. Models@Runtime overview

Algorithm Part 1

 1 DEFINITIONS Message ASK VECTORCLOCK, ASK MODEL, NOTIFICATION Type VectorClockEntry := <id: String, version ∈ N> Type Node // represents a node on the system Type Model // represents a configuration of the system Set Group:= {node: Node} Set IDS(g: Group):= {id: String | ∃ node: Node, node ∈ g & node.name = id} Set Neighbors(originator: Node, g:Group):= {node: Node | node ∈ g & originator ∈ g} Set VectorClock(originator: Node, g: Group):= {entry: VectorClockEntry | entry.id == originator.name} {entry1: VectorClockEntry |∃ node: Node, node != originator & entry1.id ∈ IDS(g) & node ∈ g} Set VectorClocks(originator: Node, g: Group):= {vectorClock: VectorClock(originator, g)} Main algorithm (see algorithm's Part 3).

3Fig. 3 .Fig. 4 .

 34 Fig. 3. Topology model of exp 1

Fig. 5 .Fig. 6 .

 56 Fig. 5. Network usage/node(in kbytes)

Fig. 8 .

 8 Fig. 7. Failure results

  The Node, score>, node ∈ Neighbors(localNode, g) && score ∈ N} nbFailure := {<node: Node, nbFail>, node ∈ Neighbors(localNode, g) && nbFail ∈ N}

	Algorithm Part 2 STATE g: Group ; changed: Boolean currentModel: Model // local version of system configuration localNode: Node // representation of local node currentVectorClock ∈ VectorClocks(localNode, g) On init(): vectorClock ← (localNode.name, 1) scores ← {Neighbors(localNode, g) × {0}} changed ← false On change (currentModel): ∀ n, n ∈ Neighbors(localNode, g) → send (n, NOTIFICATION) changed ← true Periodically do(): node ← selectPeerUsingScore() send (node,ASK VECTORCLOCK) On receive (neighbor ∈ Neighbors(localNode, g), NOTIFICATION): send (neighbor, ASK VECTORCLOCK) On receive (neighbor ∈ Neighbors(localNode,g), remoteVectorClock ∈ VectorClocks(neightbor, g)): result ← compareWithLocalVectorClock (remoteVectorClock) if result == BEFORE || result == CONCURRENTLY then send (neighbor, ASK MODEL) end if On receive (neightbor ∈ Neighbors(localNode,g),vectorClock ∈ Vectorclocks(neighbor, g), model)): result ← compareWithLocalVectorClock (targetVectorClock) if result == BEFORE then updateModel(model) mergeWithLocalVectorClock(vectorClock) else if result == CONCURRENTLY then resolveConcurrently(vectorClock, model) end if On receive (neighbor ∈ Neighbors(localNode,g), request): if request == ASK VECTORCLOCK then checkOrIncrementVectorClock() send (neighbor, currentVectorClock) end if if request == ASK MODEL then checkOrIncrementVectorClock() send (neighbor, <currentVectorClock,currentmodel>) scores := {<node: Algorithm Part 3 ALGORITHM end if

  Algorithm Part 4 SelectPeer

	Function selectPeerUsingScore()
	minScore := ∞ ; potentialPeers := {}
	for node → Neighbor(localNode, g) do
	if node ! = localNode && getScore(node) < minScore then
	minScore := getScore(node)
	end if
	end for
	for node → Neighbor(localNode, g) do
	if node ! = localNode && getScore(node) == minScore then
	potentialPeers := potentialPeers	{node}
	end if
	end for
	node := select randomly a node from potentialPeers
	updateScore(node)
	return node
	Function getScore(node ∈ Neighbors(localNode, g))
	return scores(node)
	Function updateScore(node ∈ Neighbors(localNode, g))
	oldScore := getScore(node)
	scores := scores	{node, oldScore + 2 * (nbFailure + 1)} \ {node, oldScore}

http://project-voldemort.com

More details on this project can be found in http://kevoree.org/related_projects

Acknowledgment The research leading to these results has received funding from the European Community's Seventh Framework Programme FP7/2007-2013 under grant agreement 215483 (S-Cube).


Xiaoning Guo
email: xyguo@sxicc.ac.cn

Pascal Brault
email: pascal.brault@univ-orleans.fr

Guojuan Zhi

Amaël Caillard

Guoqiang Jin

Xiangyun Guo

Structural Evolution of Plasma Sputtered Core-shell Nanoparticles for Catalytic Combustion of Methane

Keywords: Core-shell bimetallic catalysts, structural evolution, methane combustion, plasma sputtering deposition

Fe@Pd, Fe@Pt and Fe@Au core-shell nanoparticles supported by silicon carbide have been prepared by plasma sputtering deposition, and employed as the catalyst for methane combustion. The core-shell catalysts exhibit higher activities than single metallic catalysts due to surface alloying effects. With the surface alloying of the core-shell nanoparticles, Pd-O and Pt-O bonds become weak because the increasing of electron cloud density around Pd and Pt atoms due to the electron transfer from surface Fe to Pd or Pt atoms. Therefore, the activities of Fe@Pd/SiC and Fe@Pt/SiC increase with the reaction time. Whereas the activity of Fe@Au/SiC keeps invariant in the reaction due to the Fe@Au core-shell structure has high stability. Transmission electron microscopy and X-ray photoelectron spectroscopy results further confirm the structural evolution.

Introduction

Core-shell bimetallic nanoparticles have been extensively studied and widely applied in the heterogeneous catalysis, such as CO oxidation, oxygen reduction reaction and ethanol oxidation [1][2][3][4][5][6][7] . Usually, the shell layers of core-shell nanoparticles are noble metals; therefore the core-shell structures are ideal architectures for reducing dosage of noble metals and improving activity of catalysts. However, very few studies are devoted to investigating the stability of the core-shell structured nanoparticles in catalytic reactions especially under high-temperatures 8 . Moreover it is argued that catalytic core-shell nanoparticles certainly are superior to other structures such as alloys. For example, Mizukoshi et al. claimed that Au/Pd random alloy supported on TiO 2 exhibited higher activity than TiO 2 supported Au@Pd core-shell nanoparticles in photochemical evolution of H 2 from ethanol aqueous solutions under UV illumination 9 .

Various preparation techniques, such as sequential synthesis and arc discharge, have been proposed for preparation of core-shell structural catalysts [10][11][12][13] . Recently, plasma sputtering deposition has been studied extensively because the technique is feasible and no post-treatment is needed 14,15 . Considerable studies are devoted to preparation of bimetallic nanoparticles, suggesting that the plasma sputtering deposition has obvious advantages in controlling the composition of bimetallic nanoparticles [16][17][18] .

Catalytic combustion of lean natural gas and air mixture has attracted considerable attention because the process operates at low temperature and produces little NO x emission [19][20][21] . Supported noble metals (Pd, Pt and Au) have been found to have high catalytic activities for the methane combustion 19,20,22 . However noble metals are very rare and expensive. Because the activities of different catalysts are mainly related to the surface composition of metallic particles, localizing noble metals as a thin shell on common metal cores is expected to not only reduce the dosage of noble metals but also enhance the catalytic activity.

In this work, we have prepared silicon carbide supported Fe@Pd, Fe@Pt and Fe@Au core-shell bimetallic nanoparticles using the plasma sputtering deposition method and investigated the compositional and structural evolution of the different core-shell nanoparticles by transmission electron microscopy and X-ray photoelectron spectroscopy. The variation in the catalytic performance of the core-shell bimetallic catalysts reveals the structural evolution of core-shell nanoparticles in the reaction.

Experiments

Catalyst Preparation

Fe, Pd, Pt and Au sputtering depositions were performed in a cylindrical stainless steel low pressure transformer coupled plasma device. The experimental setup dedicated to the study has been described elsewhere 23,24 . Briefly, argon plasma was created in a chamber using a planar external RF antenna (13.56 MHz, 300W). The substrate used in this work is high surface area SiC 25,26 , which is prepared from a sol-gel and carbothermal reduction route and has a specific surface area of 50.8 m 2 g -1 and a pore volume of 0.13 cm 3 g -1 . The plasma depositions were performed under P = 0.5 Pa and V b = -200 V, where P is the Ar pressure (the base pressure is 5×10 -5 Pa) and V b is the bias voltage of metal targets. All the single metal targets are analytically pure (99.99%). The side length of Fe and Pt square targets is 25 mm; and the diameter of Pd and Au circular targets is 25 mm. For preparing core-shell structured nanoparticles, the deposition time was firstly set up as 20 minutes for Fe deposition and then 7.5, 7.5 and 3.75 minutes for Pd, Pt and Au deposition respectively. The corresponding samples are marked as Fe@Pd/SiC, Fe@Pt/SiC and Fe@Au/SiC. For comparison, single metallic catalysts were also prepared. The deposition time for Fe, Pd, Pt and Au is 50, 12.5, 12.5 and 6.25 minutes respectively. The corresponding samples are marked as Fe/SiC, Pd/SiC, Pt/SiC and Au/SiC.

Catalytic test

The catalytic performance of different catalysts for methane combustion was carried out in a fixed-bed quartz reactor with an inner diameter of 8 mm at atmospheric pressure, and the mixture of O 2 (20%)/CH 4 (2%)/N 2 (78%) was used as the feedstock. 100 mg of the catalyst was packed between two layers of quartz wool. The hourly space velocity was controlled to be 25000 h -1 . Since the deactivation of the SiC-supported catalysts usually demand a long time, a repeated heating-then-cooling cycle method was employed to estimate the stability of different catalysts 18,27,28 . In this method, the catalyst was programmed heated to a temperature at which the reaction obtained a near 100% methane conversion at a rate of 5 o C/min. In the heating process, the methane conversion was measured at different temperature.

Afterward, the reactor was cooled down to the temperature at which the catalyst just became inactive, and then the next reaction cycle began again. Usually, one cycle need 4 hours to finish. The composition of effluent gases was analyzed by GC-14B gas chromatograph with TDX-01 column (used to separate H2, CO, CO2, O2 and CH4) and TCD detector.

Catalyst characterization

The loadings of different metals in the catalysts were determined by Perkin-Elmer ELAN 5000 inductively coupled plasma-mass spectroscopy (ICP-MS) instrument.

The microstructures of the catalysts were analyzed by using JEM-2010 high-resolution transmission electron microscope (HRTEM). X-ray photoelectron spectroscopy (XPS) spectra of the catalysts used at different temperatures were measured on an ESCALAB 3 MKII de VG spectrometer by using Mg Kα (15kV, 20mA) X-ray source.

Results and discussion

Catalytic performances

The metal loadings of different core-shell catalysts measured by ICP-MS are listed in Table 1. In the core-shell catalysts, the Fe loadings are around 0.18 wt.% and the noble metal loadings are around 0.33 wt.%. The metal loadings of pure-metal catalysts were also determined by ICP-MS, and are around 0.5 wt% (Table 2). From Table 1 and2 Fig. 1 shows the catalytic performances of methane combustion on Fe/SiC, Pd/SiC and Fe@Pd/SiC. From Fig. 1a, Fe/SiC is active only above 650 o C; the methane conversion could reach 100% at a temperature as high as 780 o C. The activity of Pd/SiC is far better than Fe/SiC. The methane conversion on Pd/SiC increases with the reaction temperature in a relatively low temperatures range. The methane conversion on Pd/SiC begins to decline when the temperature is higher than 550 o C mainly due to decomposition of the active phase PdO into Pd 29,30 . When employing Fe@Pd/SiC as the catalyst, the methane combustion occurs at about 320 o C and the methane conversion can reach up to 100% at 530 o C. Such low loadings of Pd (0.3%) in the Fe@Pd core-shell catalyst can thus exhibit higher catalytic activity than Pd/SiC (0.5%), indicating that the core-shell structure has played an important role in the catalytic reaction. It is widely considered that the excellent catalytic performance of core-shell nanoparticles could be attributed to near surface alloy effects 31 . In the present catalyst, subsurface Pd and Pd-Fe interface layers could effectively affect the binding of adsorbates (CH 4 or O 2 ) to the Fe@Pd nanoparticle surface; and then the changes in binding enthalpies can enhance the activity of the catalyst in methane catalytic combustion process. The activity of Fe@Pd/SiC begins to decline at 570 o C which is slightly higher than that of Pd/SiC (550 o C), indicating that the Fe@Pd structure can effectively delay the declination temperature of Pd-based catalysts. Fig. 1b shows the stability of Fe@Pd/SiC catalyst in the combustion of methane.

Fe@Pd/SiC can keep the methane conversion at almost 100 % after 10 reaction cycles, indicating that high surface area SiC supported catalysts have excellent stability in methane combustion because the high thermal conductivity and chemical stability of the high surface area SiC, which can effectively hinder the sintering of active nanoparticles and thus stabilize them 18,28 . In addition, it can be obviously seen that the methane conversion of Fe@Pd/SiC at low temperatures (lower than the temperature for 100% methane conversion) increases with the cycling times. It is widely known that the activity and selectivity of catalysts extremely depend on their structures 8 . Therefore, the increase in the catalytic activity of Fe@Pd/SiC mainly results from the structural evolution of Fe@Pd core-shell nanoparticles in the catalytic process.

Fig. 2 shows the catalytic performances of Pt-and Au-based catalysts for methane combustion. Fig. 2a and 2b indicate that the core-shell catalysts exhibit slightly higher catalytic activity than their corresponding single metallic catalysts. Fig. 2c and2d show the stabilities of Fe@Pt/SiC and Fe@Au/SiC catalysts in the catalytic combustion of methane. Both catalysts can reach a 100% methane conversion and remain unchanged after 10 reaction cycles. These indicate that the high surface area SiC catalyst support can effectively stabilize metallic nanoparticles. The activity of Fe@Pt/SiC increases with the cyclic times suggesting that the Fe@Pt core-shell structure may produce similar changes as Fe@Pd/SiC during the reaction (Fig. 2c).

However, the activity of Fe@Au/SiC remains almost unaffected, indicating that the Fe@Au structure can be very stable during the reaction (Fig. 2d).

HRTEM characterization

Fig. 3 shows HRTEM images of fresh and used Fe@Pd/SiC catalysts. From Fig. 3a, the Fe@Pd core-shell structure can be clearly observed. Isolated Fe@Pd nanoparticles are highly dispersed on the SiC substrate with an average size of about 4.2 nm (Table 1). The lattice spacing of the shell layer metal is around 0.224 nm, which is indexed as Pd (111) plane. The Fe cores have an average diameter of about 2 nm and are coated with a 1.1 nm-thick Pd layer (Table 1). From Fig. 3b, the nanoparticles are still uniformly dispersed on the SiC support after reaction. However, the Fe@Pd core-shell structure has obviously changed into FePd alloy structure. The lattice spacing of 0.222 nm is indexed as the (111) plane of FePd alloy 32 . For nanoparticles which contain at least two elements, the surface atoms usually tend to lower Gibbs free energy to keep the stability of particles 33 . In our sputtering experiment, Fe atoms are firstly deposited onto the SiC surface to form nanosized Fe cores. Then Pd atoms are deposited onto the Fe-deposited SiC substrate, and then coat the Fe cores due to stronger Pd-Fe interaction than that of Pd-SiC 32 . However, the Gibbs free energy of Fe (~ 241 mJ/m 2 ) is lower than that of Pd (~ 272 mJ/m 2 ) 34 . Therefore, Fe atoms can continually move to the surface to lower the surface free energy of Fe@Pd nanoparticles, and form surface FePd alloy during the cyclic reaction. Fig. 4 shows HRTEM images of fresh and used Fe@Pt/SiC and Fe@Au/SiC catalysts. In these figures, the lattice spacings of nanoparticles are around 0.226, 0.223 and 0.234 nm, which correspond to Pt (111) plane, PtFe (111) plane and Au (111) plane, respectively 35 . The nanoparticles of different catalysts are still well dispersed. From Fig. 4a and4b, it is obvious that the Fe@Pt core-shell structure has changed into PtFe alloy after reaction because the Gibbs free energy of Fe is lower than Pt (~ 299 mJ/m 2 ) 34 . However, the Fe@Au core-shell structure can stably exist even through 10 cyclic reactions, due to the Gibbs free energy of Au (~ 177 mJ/m 2 ) being lower than Fe 34 . From Table 1 and the HRTEM images, the average sizes of metal nanoparticles in all catalysts only have a slight increase after 10 reaction cycles, indicating that high surface area SiC based catalysts have excellent stability.

XPS characterization

The XPS spectra of fresh and used Fe@Pd/SiC are shown in Fig. 5 and the corresponding binding energies (BE) of each metal and the derived surface atomic ratios are presented in Table 3. Generally, the BE values of metallic Pd are in the range 334.7-335.5 eV for Pd 3d 5/2 and 340.3-340.8 eV for Pd 3d 3/2 ; whereas the BE values of Pd 2+ are in the range 337.5-337.8 eV for Pd 3d 5/2 and 342.1-342.6 eV for Pd 3d 3/2 [36][37][38] . From Fig. 5a, The Pd BE values of fresh Fe@Pd/SiC are 335.8 and 341.2 eV, which are higher than the BE values of metallic Pd. The higher binding energies is the characteristic of Pd species with valence between 0 and +2 39 . From Fig. 5b, the concentration of Pd 0 is greater than Pd 2+ after reaction. This is due to the decomposition of PdO into Pd. However, the BE values of Pd 0 and Pd 2+ are lower than those of fresh Fe@Pd/SiC, mainly resulting form the core-shell structure being partly changed into FePd alloy. Fe atoms with lower Gibbs free energy continually migrate from core to Fe@Pd nanoparticle surface in the course of the cyclic reaction.

Meanwhile, the introduction of Fe may increase the surrounding electron cloud density around Pd 40,41 . As a result, the corresponding peaks of Pd 0 and Pd 2+ shift to lower binding energies. In addition, the main active phase of Pd-based catalysts for methane combustion is PdO, and the intensity of Pd-O bond determines the activity of catalyst 19,20 . The Pd-O bond can become weak when the electron cloud density around Pd atoms increases. Thus the CH 4 dissociation species on active sites are more easily oxidized by dissociated O species. Moreover, the products are also more easily desorbed from the active sites. As a result, the activity of Fe@Pd/SiC increases with the cyclic times.

The Fe 2p 3/2 and Fe 2p 1/2 BE values at 706.6 and 719.5 eV are attributed to metallic Fe; while the Fe 2p 3/2 and Fe 2p 1/2 at 711.3 and 723.9 eV are indexed to Fe 3+ 42,43,44 .

Therefore, the Fe component in fresh and reacted Fe@Pd/SiC exists as Fe 0 and Fe 3+ .

However, the intensity of Fe 2p peaks of used Fe@Pd/SiC is larger than the fresh one mainly because large quantity Fe atoms have migrated from cores to the surface of Fe@Pd nanoparticles. The surface atomic ratios shown in Table 3 further confirmed this result.

Fig. 6 and Fig. 7 show the XPS spectra of fresh and used Fe@Pt/SiC and Fe@Au/SiC. The corresponding BE values are also listed in Table 3. The main Pt component of used Fe@Pt/SiC is Pt 0 because PtO and/or PtO 2 are easily decomposed into Pt at higher temperatures. The Pt BE values of used Fe@Pt/SiC are lower than those of the fresh one because the Fe@Pt core-shell structure has changed into FePt alloy. The main active phases in Pt-based catalysts for methane combustion are metallic Pt, which is different with Pd-based catalysts 20 . The increase of surface Fe atoms significantly weakens the Pt-O bond 45 . Therefore the gradually weakening of Pt-O bonds results in the increasing activity of Fe@Pt/SiC with the cyclic times. By comparing with the Fig. 6c and6d, it is obvious that the intensities of Fe 2p peaks in used of reacted Fe@Pt/SiC catalyst are also larger than those of the fresh one. No obvious changes can be seen from Fig. 7, suggesting the Fe@Au core-shell structure does not change during the reaction.

Fig. 8 shows the XPS spectra of O 1s in fresh and used Fe@Pd/SiC, Fe@Pt/SiC and Fe@Au/SiC. The O 1s BE value at about 529.1 eV is attributed to the lattice oxygen associated with metal oxides, whereas the O 1s BE value at about 531.6 eV is attributed to absorbed oxygen 46 . Compare Fig. 8a and8b, it can be seen that the BE values of O1s in Fe@Pd/SiC slightly shift to lower binding energies after reation. This further confirm that the increasing of surrounding electron cloud density around Pd decline the intensity of Pd-O bond and then lower the BE values of lattice oxygen.

From Fig. 8c and 8e, only absorbed oxygen exited in fresh Fe@Pt/SiC and Fe@Au/SiC, suggesting that the metal components in both two fresh catalysts are metallic phase. In Fig. 8d, the lattice oxygen Fe@Pt/SiC in used can be attributed to Fe-O; while the BE value of adsorbed oxygen has a little decrease mainly result from the Pt-O continually weaken during the reaction. There is no change of BE values in Fe@Au/SiC has been found from the Fig. 8f, indicating that the Fe@Au core-shell nanoparticles are very stable in the methane catalytic combustion reaction.

Conclusion

Silicon carbide supported Fe@Pd, Fe@Pt and Fe@Au core-shell bimetallic catalysts were prepared by plasma sputtering deposition for methane catalytic combustion. All the core-shell catalysts exhibit higher activities than corresponding single metallic catalysts due to the near surface alloying effects. With the structural evolution of Fe@Pd and Fe@Pt nanoparticles from core-shell structure to alloy, the Pd-O and Pt-O bonds become weaker and weaker. As a result, the catalytic activities of Fe@Pd/SiC and Fe@Pt/SiC can increase gradually with the cyclic reaction times.

XPS and HRTEM results further confirm that the surface Fe concentration of Fe@Pd/SiC and Fe@Pt/SiC continually increase and the core-shell structure gradually change into alloy structure during the reaction. Whereas the activity of Fe@Au/SiC is very stable due to the stability of the Fe@Au core-shell structure in the reaction.

Tables

Table 1 The loadings of different metals, the average size of nanoparticles and the shell thickness of Fe@Pd/SiC, Fe@Pt/SiC and Fe@Au/SiC. (f) Used Fe@Au/SiC Fig. 8 XPS spectra of O 1s levels of fresh and used Fe@Pd/SiC (a and b), Fe@Pt/SiC (c and d) and Fe@Au/SiC (e and f).

  , the loadings of Pd, Pt and Fe in different catalysts have slight decreases after reaction because the volatilization of corresponding metal oxide.

Figure captions:Fig. 1

 1 Figure captions:

Fig. 2 Fig. 3 Fig. 5 XPS

 235 Fig.2 Catalytic activities of Pt/SiC and Fe@Pt/SiC (a), Au/SiC and Fe@Au/SiC (b)

Fig. 6

 6 Fig.6 XPS spectra of Pt 4f (a and b) and Fe 2p (c and d) levels of fresh and used

Fig. 7 Fig. 8 XPS

 78 Fig.7 XPS spectra of Au 4f (a and b) and Fe 2p (c and d) levels of fresh and used

Fig. 1 Fig. 3 Fig. 6 XPS

 136 Fig.1Catalytic activities (a) of Fe/SiC, Pd/SiC and Fe@Pd/SiC catalysts, and stability (b) of Fe@Pd/SiC for methane combustion.

Fig. 7

 7 Fig.7 XPS spectra of Au 4f (a and b) and Fe 2p (c and d) levels of fresh and used Fe@Au/SiC.

Table 2

 2 The loadings of different metals of Fe/SiC, Pd/SiC, Pt/SiC and Au/SiC.

		Metal loading (ICP-MS) (wt. %)		Nanoparticle	Shell
	Sample	Fe	Pd	Pt	Au	size (TEM) (nm)	thickness (nm)
		fresh used fresh used fresh used fresh used fresh used fresh used
	Fe@Pd/SiC 0.18 0.17 0.31 0.30			4.2	4.5	1.1
	Fe@Pt/SiC 0.16 0.16		0.33 0.31		3.4	3.6	0.9
	Fe@Au/SiC 0.18 0.17		0.33 0.33	4.5	4.8	1.2	1.3
			Metal loading (ICP-MS) (wt. %)
		Sample	Fe	Pd	Pt	Au
			fresh used fresh used fresh used fresh used
		Fe/SiC 0.52 0.51
		Pd/SiC		0.49 0.48
		Pt/SiC		0.50 0.47
		Au/SiC			0.49 0.49

Table 3

 3 Binding energies (eV) of Pd 3d, Pt 4f, Au 4f and Fe 2p and the surface atomic ratios of different fresh and used catalysts.

			Binding energy (eV)				Atomic ratio
	Sample	Pd	Pd	Pt	Pt	Au	Au	Fe	Fe	Pd (Pt or Au)
		3d 5/2	3d 3/2	4f 7/2	4f 5/2	4f 7/2	4f 5/2	2p 3/2	2p 1/2	/ Fe
	Fresh Fe@Pd/SiC	335.8 337.6 341.2 342.6					706.5 719.5	69:31
	Used Fe@Pd/SiC	335.4 337.0 340.7 342.3					711.3 723.9	51:49
	Fresh Fe@Pt/SiC			71.6 75.2			706.6 719.5	70:30
	Used Fe@Pt/SiC			71.1 74.8			711.4 723.8	46:54
	Fresh Fe@Au/SiC					84.1 87.8 706.6 719.7	73:27
	Used Fe@Au/SiC					84.0 87.9 711.3 724.0	71:29

Acknowledgments

The work was financially supported by NSFC (Ref. 20973190), the in-house research project of SKLCC (Ref. SKLCC-2008BWZ010), and the National Basic