E D Andjel

CONVERGENCE IN DISTRIBUTION FOR SUBCRITICAL 2D ORIENTED PERCOLATION SEEN FROM ITS RIGHTMOST POINT

Keywords: Oriented percolation, rightmost point, subcritical

We study subcritical two-dimensional oriented percolation seen from its rightmost point on the set of infinite configurations which are bounded above. This a Feller process whose state space is not compact and has no invariant measures. We prove that it converges in distribution to a measure which charges only finite configurations.

Introduction and main results.

1.1. Background. Two-dimensional oriented percolation and its continuous time analog the one-dimensional contact process, seen from their rightmost point, have been studied in several papers. Durrett [START_REF] Durrett | Oriented percolation in two dimensions[END_REF] proved that in the critical and supercritical phase there exists an invariant measure. Then, Schonmann proved that there are no such measures in the subcritical phase [START_REF] Schonmann | Absence of a stationary distribution for the edge process of subcritical oriented percolation in two dimensions[END_REF]. These two papers consider only the discrete time model, but their results hold also for some continuous time one-dimensional process which include the contact process (see [START_REF] Andjel | Edge processes of one-dimensional stochastic growth models[END_REF]). Galves and Presutti [START_REF] Galves | Edge fluctuations for the one-dimensional supercritical contact process[END_REF] proved that the onedimensional contact process seen from the rightmost point converges in the supercritical phase to a unique invariant measure. This last result was then extended by Cox, Durrett and Schinazi [START_REF] Cox | The critical contact process seen from the right edge[END_REF] to the critical phase. There are no difficulties in adapting these convergence results to the discrete time setting. Finally, we mention [START_REF] Galves | Edge fluctuations for the one-dimensional supercritical contact process[END_REF] and [START_REF] Kuczek | The central limit theorem for the right edge of supercritical oriented percolation[END_REF] where the position of the rightmost point is shown to satisfy a central limit theorem. In this paper, we prove that convergence of the discrete time process seen from the rightmost point also occurs in the subcritical phase although there are no invariant measures. 1.2. Definitions. Let Λ = {(x, y) : x, y ∈ Z, y ≥ 0, x + y ∈ 2Z}.

(1.1)

E. D. ANDJEL

Draw oriented bonds from each point (m, n) in Λ to (m + 1, n + 1) and to (m -1, n + 1). In this paper, we suppose that bonds are open independently of each other and that each bond is open with probability p ∈ (0, 1). To formalise this, let B be the set of all bonds with both endpoints in Λ and assume that {Z b : b ∈ B} is a collection of i.i.d. random variables whose distribution is a Bernoulli with parameter p. A bond b will be considered open (closed) if Z b = 1 (Z b = 0). The event consisting on the existence of an open path from A to B, where A and B are subsets of Λ, will be denoted by {A → B} and its complement by {A B}. When either A or B or both are singletons, say {x} and {y}, respectively, we will write {x → y}, {x → B}, etc.

Given a subset A of 2Z we let

ξ A n = {y : (y, n) ∈ Λ and (x, 0) → (y, n) for some x ∈ A}, (1.2) n = 0, 1, . . . .
Then, (ξ A n , n ≥ 0) is a Markov chain taking values in the subsets of 2Z at even times and of 2Z + 1 at odd times.

Let A be an infinite subset of 2Z such that sup A < ∞. Then, for all n > 0, the supremum of ξ A n is finite and a simple Borel-Cantelli argument shows that ξ A n is a.s. infinite. For such initial conditions, we let

r(ξ A n ) = sup ξ A n and ζ A n = {x -r(ξ A n ) : x ∈ ξ A n }. Then (ζ A n , n ≥ 0
) is a Markov chain on infinite subsets of 2Z -=: {0, -2, -4, . . .} containing 0. For finite subsets A we may also define the Markov chain (ζ A n , n ≥ 0) by simply adopting the convention:

ζ A n = ∅ if ξ A n = ∅.
Obviously, ∅ is an absorbing state for both (ξ A n , n ≥ 0) and (ζ A n , n ≥ 0). In the sequel,

S = {infinite subsets of 2Z -containing 0}, (1.3) S 0 = {finite subsets of 2Z -containing 0} (1.4) and S = {subsets of 2Z -}. (1.5)
We will consider S and S 0 as subsets of S which we identify with {0, 1} 2Z - by means of the bijection: F (A) = 1 A . Then, S inherits the product topology of {0, 1} 2Z -and becomes a compact space. The subsets S and S 0 of S are now endowed with the induced topology. Probability measures on either S or S 0 will be seen as measures on S and the space of all probability measures on S will be endowed with the topology of weak convergence.

Standard coupling arguments show that P (ξ 0 n = ∅ for all n) increases with p and we can define the critical value p c of the parameter p as the supremum of its values for which the above probability is 0. It is well known (see [START_REF] Durrett | Oriented percolation in two dimensions[END_REF]) that 0 < p c < 1. Throughout this paper, we assume that p ∈ (0, p c ). 1.3. Theorems. Before stating our results, we recall that a quasi-stationary distribution of a Markov chain (X n ; n ≥ 0) on S 0 ∪ {∅} with absorbing state ∅ is a probability measure ν on S 0 such that P ν (X n = x|T > n) = ν(x) for all n ∈ N and x ∈ S 0 , where T = inf{k : X k = ∅}. We refer the reader to [START_REF] Ferrari | R-positivity, quasi-stationary distributions and ratio limit theorems for a class of probabilistic automata[END_REF] for more information concerning quasi-stationary distributions. Our first theorem is not new, it is immediately obtained from Theorem 1 of [START_REF] Ferrari | R-positivity, quasi-stationary distributions and ratio limit theorems for a class of probabilistic automata[END_REF]. n given {T > n} converges as n goes to infinity to a probability measure ν on S 0 . Moreover, ν is the minimal quasi-stationary distribution of the ζ n process on S 0 ∪ {∅}.

We now state our main result which was conjectured by Galves, Keane and Meilijson. As expected by the authors of [START_REF] Ferrari | R-positivity, quasi-stationary distributions and ratio limit theorems for a class of probabilistic automata[END_REF] (see Remark 7 in page 606 of that reference), Theorem 1.1 is the key ingredient to prove it. Theorem 1.2. Suppose 0 < p < p c . Then, for any A ∈ S the distribution of ζ A n converges as n goes to infinity to ν, where ν is as in Theorem 1.1.

The paper is organised as follows: Section 2 starts explaining the strategy we will follow, continues stating two lemmas and then deduces Theorem 1.2 from these lemmas. Then, in Section 3 we prove those two lemmas.

2. Proof of Theorem 1.2. We start this section introducing some more notation: Let f be real-valued function defined on S ∪ S 0 . We say that f is a cylinder function depending only on coordinates -2r, . . . , -2 if there exists a function g defined on subsets of {-2r, . . . , -2} such that f (A) = g(A ∩ {-2r, . . . , -2}) for all A ∈ S ∪ S 0 . For such functions, will denote the supremum norm

f = sup A∈S∪S 0 |f (A)|. For (x, m) ∈ Λ let C x,m = {(y, k) ∈ Λ : k ≥ m, |y -x| ≤ k -m}
and call this set the cone emerging from (x, m). For r ∈ N, call level r the set

L r = {(x, n) ∈ Λ : n = r}.
We will say that level n is higher than level m if n ≥ m. In the sequel, ν r will be the distribution of ζ 0 r given {T > r} where T = inf{k : ζ 0 k = ∅} and A will be fixed but arbitrary element in S.

E. D. ANDJEL

We now sketch the proof of Theorem 1.2: We first find the rightmost point x 0 of A satisfying ξ x 0 n = ∅. We would like to apply Theorem 1.1 but cannot do it immediately because we are conditioning not only on {ξ x 0 n = ∅} but also on {ξ y n = ∅} for all y ∈ A ∩ {z : z > x 0 }. However since p < p c , there is a positive probability that no point of C x 0 ,0 can be attained from {(y, 0) : y > x, y ∈ A} following open paths. If this occurs, then the distribution of ζ x 0 n is ν n . If this fails to happen we look at the highest level in C x 0 ,0 attained from {(y, 0) : y > x 0 , y ∈ A} and repeat the argument from that level. We keep doing so until the corresponding emerging cone is not attained. Once this happens, we will derive from p < p c that the elements of ξ A n \ ξ x 0 n are far to the left of ξ x 0 n for large n. In carrying out this approach, the main difficulty comes from keeping track of several conditionings. To make this argument rigorous, we begin defining two sequences of r.v.'s Y i and X i as follows: Let

Y 0 = 0 (2.1)
and

X 0 = sup{x ∈ A : (x, 0) → L n }. (2.2)
Then, given Y 0 , Y 1 , . . . , Y i and X 0 , . . . , X i , we let

Y i+1 = sup{k : ∃u, v : (u, Y i ) → (v, k) (2.3) with u > X i , u ∈ ξ A Y i and (v, k) ∈ C X i ,Y i } ∨ Y i and X i+1 = sup{x : x ∈ ξ A Y i+1 and (x, Y i+1 ) → L n }. (2.4)
Note that Y i ≤ Y i+1 and that as soon as Y i = Y i-1 both sequences become constant. The reader may find helpful to have now a first look at Figure 1 in the next section.

We now state two lemmas which will be proved in the next section. In the second of these lemmas, we use the fact that on the event {0 → L n } the process (ζ 0 k , k = 0, . . . , n) takes values on S 0 .

Lemma 2.1.

Let I = inf{i : Y i = Y i+1 }.
Then, there exists β > 0 such that for all m:

(a) P (I ≥ m) ≤ exp(-βm), (b) P (Y I ≥ m 2 ) ≤ (m + 1) exp(-βm) and (c) P (I ≤ m, Y I ≤ m 2 ) ≥ 1 -(m + 2) exp(-βm).
Lemma 2.2. Let A be an element of S, let f be a cylinder function on S ∪ S 0 depending only on the coordinates -2r, . . . , -2 and let β be as in Lemma 2.1. Then, for all i, j ≤ n,

|E(f (ζ A n )|I = i, Y i = j) -E(f (ζ 0 n )|0 → L n )| ≤ 2 f (n + r) exp(-β(n -j)).
We now proceed to prove our main result.

Proof of Theorem 1.2. Let f be a cylinder function on S ∪ S 0 depending only on the coordinates -2r, . . . , -2 and let m = ⌊n 1/3 ⌋ where ⌊•⌋ denotes the integer part. By part (c) of Lemma 2.1, we have

E(f (ζ A n )) - m i=0 m 2 j=0 E(f (ζ A n )|I = i, Y i = j)P (I = i, Y i = j) ≤ f (m + 2) exp(-βm).
Therefore,

|Ef (ζ A n ) -E(f (ζ 0 n )|0 → L n )| ≤ m i=0 m 2 j=0 (|E(f (ζ A n )|I = i, Y i = j) -E(f (ζ 0 n )|0 → L n )|)P (I = i, Y i = j) + 2 f (1 -P (I ≤ m, Y I ≤ m 2 )) ≤ m i=0 m 2 j=0 (2 f (n + r) exp(-β(n -j))P (I = i, Y i = j)) + 2 f (1 -P (I ≤ m, Y I ≤ m 2 )) ≤ 2 f (n + r) exp(-β(n -m 2 )) + 2 f (1 -P (I ≤ m, Y I ≤ m 2 )),
where the second inequality follows from Lemma 2.2. Since m = ⌊n 1/3 ⌋ this and part (c) of Lemma 2.1, imply that

lim n |Ef (ζ A n ) -E(f (ζ 0 n )|0 → L n )| = 0,
and the result follows from Theorem 1.1.

3. Proofs of Lemmas 2.1 and 2.2. In this section, for r ∈ N and x ∈ 2Z -, G r

x will denote the σ-algebra generated by the random variables which determine the state of the bonds with both vertices in (

-x/2 i=0 C x+2i,0 ) ∩ ( r j=0 L j )
, G r will denote the σ-algebra generated by the random variables which determine the state of the bonds with both vertices in r i=0 L i and G ′ r will denote the σ-algebra generated by the random variables which determine the state of the bonds with both vertices in ∞ i=r L i . Besides this, an event belonging to a σ-algebra generated by random variables determining the state of a finite number of bonds will be called an elementary cylinder of that σ-algebra if it is nonempty and does not contain any nonempty proper subset of that σ-algebra. The first lemma of this section is an immediate consequence of the exponential decay of P (ξ x n = ∅) (see Section 7 of [START_REF] Durrett | Oriented percolation in two dimensions[END_REF]) and we omit its proof. In our next lemma, for notational convenience we let r 0 = 0 and recalling (2.1)-(2.4), consider events of the form

G(x 0 ) = {X 0 = x 0 } and for i ≥ 1 G(r 1 , . . . , r i ; x 0 , . . . , x i ) = {Y 1 = r 1 , . . . , Y i = r i , X 0 = x 0 , . . . , X i = x i },
where 0 ≤ r 1 ≤ • • • ≤ r i are integers and (x 0 , 0), (x 1 , r 1 ), . . . , (x i , r i ) ∈ Λ.

Lemma 3.2. Let i be a nonnegative integer and let F be an elementary cylinder in G r i

x 0 having a nonempty intersection with G(r 1 , . . . , r i ; x 0 , . . . , x i ). Then, L r i contains i + 2 finite subsets A i,1 , . . . , A i,i , B i , D i determined by F, n, x 0 , . . . , x i , r 1 , . . . , r i only and such that

F ∩ G(r 1 , . . . , r i ; x 0 , . . . , x i ) = F ∩ {(x i , r i ) → L n } ∩ {u L n ∀u ∈ D i } (3.1) ∩ {u L r i +1 ∀u ∈ B i } ∩ i j=1 {u L n , u C x j-1 ,r j-1 ∀u ∈ A i,j } . Moreover, ξ A r i (x i ) = 1, {z : z > x i , ξ A r i (z) = 1} = D i ∪ B i ∪ i j=1 A i,j (3.2) and x i < d i < b i < a i,i < a i,i-1 < • • • < a i,1 (3.3) ∀d i ∈ D i , b i ∈ B i , a i,j ∈ A i,j , j = 1, . . . , i.
Remark. If F is disjoint of G(r 1 , . . . , r i ; x 0 , . . . , x i ), we may extend the definition of the sets A i,1 , . . . , A i,i , B i , D i by letting them be the empty set. In this way, they become random G r i

x 0 -measurable sets, hence independent of the σ-algebra G ′ r i .

Proof of Lemma 3.2. To follow this proof, we recommend the reader to look at Figure 1. This may help visualizing the different sets involved in it. We proceed by induction on i. If i = 0, then G r i

x 0 = G 0 x 0 is the trivial σ-algebra. Hence, F must be the whole probability space and the statement holds with D 0 = {(z, 0) : z > x 0 , z ∈ A} and B 0 = ∅. Assume the statement holds for some given i, and let F ′ be an elementary cylinder of G r i+1

x 0 . Call F the unique elementary cylinder of G r i

x 0 which contains F ′ . Then, by the inductive hypothesis there are i + 2 subsets A i,1 , . . . , A i,i , B i , D i for which (3.1), (3.2) and (3.3) hold. Now, define the following subsets of L r i+1 :

A i+1,j = {(x, r i+1 ) : A i,j → (x, r i+1 )} (i ≥ 1, j = 1, . . . , i), A i+1,i+1 = {(x, r i+1 ) / ∈ C x i ,r i : D i → (x, r i+1 )} (i ≥ 0), B i+1 = {(x, r i+1 ) ∈ C x i ,r i : D i → (x, r i+1 )} and
D i+1 = {(x, r i+1 ) ∈ C x i ,r i : x > x i+1 , (x i , r i ) → (x, r i+1 )} \ B i+1 .
It is now tedious but straightforward to verify that these sets satisfy (3.1), (3.2) and (3.3) with F ′ and i + 1 replacing F and i, respectively. Proposition 3.1. Let β be as in Lemma 3.1. Then, for all m, i and all x 0 , . . . , x i , r 1 , . . . , r i we have:

P (Y i+1 -Y i ≥ m|G(r 1 , . . . , r i , x 0 , . . . , x i )) ≤ exp(-βm). (3.4) E. D. ANDJEL
Proof. Call Υ the set of all paths from (x i , r i ) to L n . Given a path γ ∈ Υ, call A γ,r the oriented graph composed by the bonds having both vertices in n j=r i L j and at least one vertex strictly to the right of γ, and by the vertices of such bonds. Let Γ r be the rightmost open path starting from (x i , r i ) and attaining L n . Let γ be a possible value of Γ r . Note that the event {Γ r = γ} is constituted by the configurations for which γ is open but there is no open path from a vertex of γ contained in A γ,r and reaching either L n or another point in γ. Hence, the event {Γ r = γ} is the intersection of the event {γ is open} and a decreasing event D(γ) on the graph A γ,r .

Let F be an elementary cylinder in G r i x 0 having a nonempty intersection with G(r 1 , . . . , r i ; x 0 , . . . , x i ). To prove the proposition it suffices to show that for some β > 0 which depends only on p, we have

P (Y i+1 -Y i ≥ m|F ∩ G(x 0 , . . . , x i , r 1 , . . . , r i )) ≤ exp(-βm). (3.5) By Lemma 3.2 on the event F ∩ G(x 0 , . . . , x i , r 1 , . . . , r i ) ∩ {Y i+1 -Y i ≥ m}, there is a point in u ∈ D i such that u → ( n j=r i +m L j ) ∩ C x i ,r i . Hence, (3.5

) will follow if we prove

P u → C x i ,r i ∩ n j=r i +m L j for some u ∈ D i (3.6)
F ∩ G(x 0 , . . . , x i , r 1 , . . . , r i ) ≤ exp(-βm). By Lemma 3.2 this can be written as

P u → C x i ,r i ∩ n j=r i +m L j for some u ∈ D i F ∩ {(x i , r i ) → L n } ∩ {u L n ∀u ∈ D i } ∩ {u L r i +1 ∀u ∈ B i } ∩ i j=1 {u L n , u C x j-1 ,r j-1 ∀u ∈ A i,j } ≤ exp(-βm).
Since the state of the bonds above L r i is independent of F this is equivalent to

P u → C x i ,r i ∩ n j=r i +m L j for some u ∈ D i {(x i , r i ) → L n } ∩ {u L n ∀u ∈ D i } ∩ {u L r i +1 ∀u ∈ B i } ∩ i j=1 {u L n , u C x j-1 ,r j-1 ∀u ∈ A i,j } ≤ exp(-βm).
Since {(x i , r i ) → L n } is a disjoint union of the events {Γ r = γ} where γ ranges over Υ, it suffices to show that for any γ ∈ Υ we have

P u → C x i ,r i ∩ n j=r i +m L j for some u ∈ D i {Γ r = γ} ∩ {u L n ∀u ∈ D i } ∩ {u L r i +1 , ∀u ∈ B i } (3.7) ∩ i j=1 {u L n , u C x j-1 ,r j-1 ∀u ∈ A i,j } ≤ exp(-βm).
But, as explained at the beginning of this proof, the left-hand side above can be written as

P u → C x i ,r i ∩ n j=r i +m L j for some u ∈ D i {γ is open} ∩ D(γ) ∩ {u L n ∀u ∈ D i } ∩ {u L r i +1 ∀u ∈ B i } (3.8) ∩ i j=1 {u L n , u C x j-1 ,r j-1 ∀u ∈ A i,j } .
Now, let V (γ) be the set of vertices of γ. Then, noting that

u → C x i ,r i ∩ n j=r i +m L j for some u ∈ D i ∩ {u L n ∀u ∈ D i } ∩ {γ is open} (3.9) = u → C x i ,r i ∩ n j=r i +m L j within A γ,r for some u ∈ D i ∩ {u L n ∀u ∈ D i } ∩ {γ is open}, and that {u L n ∀u ∈ D i } ∩ {γ is open} (3.10) = {u L n ∪ V (γ) ∀u ∈ D i } ∩ {γ is open}, 10 E. D. ANDJEL
(3.8) can be written as

P u → C x i ,r i ∩ n j=r i +m L j within A γ,r for some u ∈ D i {γ is open} ∩ D(γ) ∩ {u L n ∪ V (γ) ∀u ∈ D i } (3.11) ∩ {u L r i +1 ∀u ∈ B i } ∩ i j=1 {u L n , u C x j-1 ,r j-1 ∀u ∈ A i,j } .
Since {γ is open} is independent of all the other events involved in the above expression, (3.11) is equal to

P u → C x i ,r i ∩ n j=r i +m L j within A γ,r for some u ∈ D i D(γ) ∩ {u L n ∪ V (γ) ∀u ∈ D i } ∩ {u L r i +1 ∀u ∈ B i } (3.12) ∩ i j=1 {u L n , u C x j-1 ,r j-1 ∀u ∈ A i,j } . Since {u → C x i ,r i ∩ ( n j=r i +m L j ) within A γ,r for some u ∈ D i } is an increas- ing event while D(γ) ∩ {u L n ∪ V (γ) ∀u ∈ D i } ∩ {u L r i +1 ∀u ∈ B i } ∩ ( i j=1 {u L n , u C x j-1 ,r j-1 ∀u ∈ A i }) is decreasing, (3.12) is bounded above by P u → C x i ,r i ∩ n j=r i +m L j within A γ,r for some u ∈ D i ≤ P u → C x i ,r i ∩ n j=r i +m L j for some u ∈ L r i (3.13)
to the right of (x i , r i ) .

The proposition now follows from Lemma 3.1.

(ζ A ′ n )|{0 → L n }) -E(f (ζ 0 n )|{0 → L n })| ≤ 2(n + r) f exp(-βn) ∀n ∈ N.
Proof. Let Φ be the (random) set of points in levels 1, 2, . . . , n which can be reached from (0, 0) following an open path. The event {0 → L n } is a disjoint union of events of the form {Φ = κ} where κ ranges over all values of Φ containing at least one point in L n . Then write Before starting the proof of Lemma 2.2, we need to introduce some further notation: T will the map sending subsets of 2Z + k into subests of 2Z + k + 1 given by T (A) = {x -1; x ∈ A} and for (x, n) ∈ Λ let G + x,n be the σ-algebra generated by the random variables determining the state of the bonds whose vertices are in n i=0 L i and by the bonds having at least one vertex strictly to the right of {(x + i, n + i); i = 0, 1, . . .} and let G - x,n be the σ-algebra generated by the random variables determining the state of all the other bonds. If B is an infinite subset of 2Z + k which is bounded above, we define for n ≥ k: ξ

|E(f (ζ A ′ n )|Φ = κ) -E(f (ζ 0 n )|Φ = κ)| ≤ 2 f

Theorem 1 . 1 .

 11 Suppose 0 < p < p c and let T = inf{n : ζ 0 n = ∅}. Then the conditional distribution of ζ 0

Lemma 3 . 1 .

 31 There exists a constant β > 0 such that for all x ∈ 2Z and all m ∈ N we haveP (y, 0) → C x,0 ∩ ∞ j=mL j for some y > x ≤ exp(-βm).

Fig. 1 .

 1 Fig. 1. Here I = 2, (X0, Y0) = (x0, 0), (X1, Y1) = (x1, r1), (X2, Y2) = (x2, r2), the doted lines are the emerging cones from these three points and the full lines are the open paths starting from A.

Proof of Lemma 2 . 1 . 1 j=0P

 211 If I ≥ m then the sequence Y 0 , . . . , Y m is strictly increasing. But it follows from Proposition 3.1 that P (Y i+1 > Y i |Y 1 , . . . , Y i ) ≤ exp(-β) a.s. Hence (a) follows by induction in i. To prove (b) write P (Y I ≥ m 2 ) ≤ P (I > m) + P (Y m ≥ m 2 ) ≤ exp(-βm) + m-(Y j+1 -Y j ≥ m) ≤ exp(-βm) + m exp(-βm), where the second inequality follows from part (a) and the last one from Proposition 3.1. Part (c) follows easily from parts (a) and (b). Lemma 3.3. Let A ′ ∈ S ∪ S 0 and let f be a cylinder function on S ∪ S 0 depending only on the coordinates -2r, . . . , -2. Then, |E(f

P

  (-2i → L n off κ) ≤ 2(n + r) f exp(-βn)and the lemma follows.

  B k,n = {z : (x, k) → (z, n) for some x ∈ B}, r(ξ B k,n ) = sup(ξ B k,n ) and ζ B k,n = {x -r(ξ B k,n ) : x ∈ ξ B k,n }. As before (ζ B k,n , n ≥ k) is a Markov chain on infinite subsets of 2Z -containing 0. Proof of Lemma 2.2. Since E(f (ζ A n )|I = i, Y i = j) is a convex combination of E(f (ζ A n )|I = i, Y i = j, X i = x i )where x i runs over all possible values of X i , it suffices to show that for allx i |E(f (ζ A n )|I = i, Y i = j, X i = x i ) -E(f (ζ 0 n )|0 → L n )| (3.14)≤ 2 f (n + r) exp(-β(n -j)).

Acknowledgements. The author wishes to thank A. Galves who mentioned to him the plausibility of Theorem 1.2, P. Ferrari for conversations concerning the results of [4], F. Ezanno who provided Figure 1 and careful reading which improved this paper, and the referees for their useful comments.

E. D. ANDJEL

But on the event {I = i, Y i = j, X i = x i } it happens that x i is the rightmost point of ξ A j from which there is an open path to L n . Therefore, (3.14) will follow if we show that for all infinite subset A ′ of L j such that sup A ′ = x i we have

.


Maher Barkallah

Slim Bouaziz
email: 2slim.bouaziz1@gmail.com

Najib Belhadj Messaoud
email: najibebm@yahoo.fr

Jean-Yves Choley
email: jean-yves.choley@supmeca.fr

Mohamed Haddar
email: mohamed.haddar@enis.rnu.tn

Vibration response simulation of a machine tool spindle system in non-stationary regime

Keywords: tool, stiffness, cutting force, spindle, fluctuation, torque

Bearing and gear are one of the most important mechanical sources for vibration and noise generation in machine tool spindles. In this paper, we study the non-linear dynamic behavior of a machine tool spindle system in transient regime. Driving and driver rotors are, respectively, powered by a motor torque and loaded by the cutting force. They are supported by two identical Rolling Bearings (RBs). Gear excitation is induced by the motor torque and load variation in addition to the fluctuation of meshing stiffness due to the variation of input rotational speed. The dynamic parameters of RBs are modeled by stiffness and damping matrices computed by the derivation of the bearing forces. The equations of motion are solved iteratively using Newmark time integration method. The numerical results of the dynamic responses of the system come to confirm the significant effect of the transient regime on the dynamic behavior of a gear set.

Introduction

Bearings play an important role in machine tool spindle systems. Compared with hydrostatic, aerostatic or magnetic bearings [START_REF] Inayat-Hussain | Geometric coupling effect on the bifurcations of a flexible rotor response in active magnetic bearings[END_REF], rolling element bearings are still most commonly used today in the spindles, which can provide the required precision, load carrying capacity, and spindle speeds. The dynamic behavior from the transient state chaotic motion to the steady state periodic and quasi periodic motions is also established by Zhang et al. [START_REF] Wei | Transient and steady nonlinear responses for a rotor-active magnetic bearings system with time-varying stiffness[END_REF]. Luisa et al. [START_REF] Villa Luisa | Statistical fault diagnosis based on vibration analysis for gear test-bench under non-stationary conditions of speed and load[END_REF] reveal how to take advantage of the information on vibrations from the mechanical system in a varied range of speed and load conditions. Acceleration signals from the accelerometers are treated using a new formulation to study the gear motions. In the developed concept of supports for aviation purposes, some simulations have been performed. A necessity of increasing bearings load capacity and damping has occurred by Kozanecka et al. [START_REF] Kozanecka | Active magnetic damper in a power transmission system[END_REF].

In rotating machine, to provide speed and torque conversions from a rotating power source to connected mechanical devices, gear reducers are often used but its represent a principal source of vibration and noise due to excitations associated with the conditions of contact between gear pair teeth. Kang et al. [START_REF] Ma Ru | Measurement of vibratory motions of gears supported by compliant shafts[END_REF] made an experimental results specify that the proposed accelerometer-based measurement methodology is operational in measuring not only the torsional motions but also the other types of translational, rotational and axial motions of a gear. To study the dynamic behavior of a spur gear pair it is required to define the stiffness in the contact zone named gear mesh stiffness. The gear set have internal and external sources of excitations. The internal one is made by the time varying mesh stiffness. These fluctuations are considered the main source of system excitation and the origin of the observed noise and vibrations [START_REF] Chaari | Analytical modelling of spur gear tooth crack and influence on gearmesh stiffness[END_REF]. Bartelmus [START_REF] Bartelmus | Mathematical modelling of gearbox vibration for fault diagnosis[END_REF][START_REF] Bartelmus | Mathematical modelling and computer simulations as an aid to gear box diagnostics[END_REF] presented the varying mesh stiffness in dynamic model of a spur gear system to investigate the different responses in presence of defects. Chaari et al. [START_REF] Chaari | Effect of spalling or tooth breakage on gear mesh stiffness and dynamic response of a one-stage spur gear transmission[END_REF][START_REF] Chaari | Influence of manufacturing errors on the dynamic behavior of planetary gears[END_REF] introduced the varying mesh stiffness in a dynamic model of a spur gear system in order to study its dynamic responses in presence of defects. Added frequency components were observed in response spectra. Sachidananda et al. [START_REF] Sachidananda | Experimental investigation of fatigue behavior of spur gear in altered tooth-sum gearing[END_REF] presented analytical and experimental methods to investigate and compare the altered tooth sum gearing against the standard tooth-sum gearing. The experiments were performed using a power recirculating type test rig. The tooth loads for the experimental investigations were determined considering the surface durability of gears.

Bartelmus [START_REF] Bartelmus | Mathematical modelling of gearbox vibration for fault diagnosis[END_REF][START_REF] Bartelmus | Mathematical modelling and computer simulations as an aid to gear box diagnostics[END_REF] and Walha et al. [START_REF] Walha | Nonlinear dynamics of a two-stage gear system with mesh stiffness fluctuation, bearing flexibility and backlash[END_REF] studied the effect of backlash on the dynamic behavior of a two stage gearbox. They concluded that tooth separation occurs at the transient regime with increase in the vibration level. Sika and Velex [START_REF] Sika | Analytical and numerical analysis of gears in the presence of engine acyclism[END_REF][START_REF] Sika | Analytical and numerical analysis of gears in the presence of engine acyclism[END_REF] proposed a model to study, in transient regime, the influence of gear tooth geometry, the backlash and the backstrike effect on the dynamic behavior of gear system. A gear set powered by an electric motor in the start time was also investigated by Hugues [START_REF] Hugues | Contribution to the Dynamic Study, During Startup, of Drivelines Gear Driven By Electric Motors[END_REF]. He noticed that the existence of distributed faults induces an increase in the vibration levels when load increases.

Model of a single stage spur gear reducer

A single stage spur gear model with eight degrees of freedom is presented in Figs. 1 and2 The system is sought by motor torque and load torque . The gear mesh stiffness ( ) and damping ( ) are modelled by linear spring and damping acting on the line of meshing teeth action (Fig. 3). The displacement along the line of action is expressed by [START_REF] Velex | Dynamic tooth loads and quasi-static transmission errors in helical gearsapproximate dynamic factor formulae[END_REF]:

( ) = ( -)sin + ( -)cos + + , (1)
where and are the translation of the pair of gear ( = 1, 2). is the angular displacement of the component in part ( , = 1, 2), is the pressure angle, and are, respectively, the base radius of the pinion and the wheel, is the rotational speed of the pinion. The mesh period can be defined by:

= 60 . (2)
The vector of the degrees of freedom can be expressed by:

= [ ] . (3)

Derivation of the equations of motion

Rolling bearing modeling

The forces and moments, exerted by the th rolling bearing supporting the spindle (Fig. 4), are defined, and using the Hertz contact theory, as:

= -Δ , cos cos sin cos , (4)
where, and represents respectively the ball angular position and the loaded contact angle, presents the radial position of the outer race curvature centre of the th ball, represents the Hertz contact constant, deduced from curvature radius of the elements in contact, and Δ represents the elastic deformations of the th ball.

This last is defined as follows:

Δ = Δ +Δ -, (5)
where, Δ and Δ present respectively the radial and axial deflections. They are expressed in the reference , , , functioning of the elastic and rigid movements of the th node of the spindle: The new loaded contact angle can be deduced:

tan = Δ Δ . ( 7
)
The scaling default is modeled by a half sine wave introduced at the th ball deflections expressions, which will be:

Δ = Δ +Δ --. (8)
In the present study the classic linearised model with four stiffness and damping coefficients is used for the RBs modeling. In this model presented in Fig. 5, the forces at each bearing are assumed to obey the governing equations of the following form [START_REF] Slim | Haddar Mohamed Dynamic behaviour of hydrodynamic journal bearings in presence of rotor spatial angular misalignment[END_REF]:

+ = , (9)
where and are the RBs degrees of freedom, [ ] and [ ] are respectively the RBs stiffness and damping matrices. They are expressed by:

= , = . (10)
The derived equations of motion can be represented in a matrix-vector form. It can be written as:

[ ] + [ ( )] + [ ( )] = ( ) , (11)
where, represents the mass matrix given by:

[ ] = 0 0 . (12)
[ ( )] includes the RBs stiffness ( , = , ), the shaft torsional stiffness , and the time varying gear mesh stiffness ( ). It is expressed by:

( ) = 0 + (
[ ( )] includes the RBs damping ( , = , ), the shaft torsional damping , and the time varying gear mesh damping ( ). It is expressed by:

Electric motor modeling

The rotational velocity of an electric motor can be related to its torque and the receiver torque by:

Ω = - , (18)
where is the equivalent moment of inertia of the rotating parts. The relation between driving torque and the motor speed have the following form [START_REF] Wright | Class Notes on Design and Analysis of Machine Elements[END_REF]:

= 1+( -) - , (19)
where is the slip, is the torque at breakdown, and are two constants characteristic of the motor and is the proportional drop in speed given by:

=1- Ω Ω , (20)
where Ω and Ω are, respectively, the synchronous speed and the actual rotational speed of the motor.

In the present study a pump is chosen as a load. Load torque characteristic is proportional to the square of the rotational speed [START_REF] Khabou | Study of a spur gear dynamic behavior in transient regime[END_REF]:

=ℎ Ω , ( 21
)
where ℎ is a coefficient depending on the driven system kind. Eq. ( 17) can be written as follows:

Ω = 1 (Ω, ) - (Ω, ) . (22)

Numerical simulations

Using Newmark integration method, the dynamic behavior of spur gear system is studied by a numerical simulation. The mechanical system is powered by an electric motor. The principal characteristic parameters of the electric motor are given in Table 1. The principal parameters of the spur gear transmission are given in Table 2. Tables 3 presents the eigenfrequencies of the system. The mechanical transmission is loaded by a machine tool having a n i n e r t i a m o m e n t = 45×10 -4 Kg•m 2 . Using the Euler algorithm to resolve Eq. ( 11), the evolution of motor rotational speed versus time is presented in Fig. 6. Two different regimes are determinate. The first is the transient regime observed in the time interval [0, ], characterized by an increasing rotational speed. The second regime is the steady state observed in a time greater than and subsequently the motor reaches its nominal rotational speed Ω . The transient regime duration is = 0.185 s. The evolution of motor and receiver torques of the studying system is shown in Fig. 7. The evolution of the meshing stiffness in transient regime is shown in Fig. 8; numerical results show that it is modeled as a non-periodic function. It should be noted that the meshing period decreases when the speed increases.

The dynamic response registered on the pinion is presented in Fig. 9(a). Two different behaviors are determinate: The first period ( ) corresponding to the interval of time [0, ] in which the displacement is variable and non-periodic with high vibration amplitudes (Fig. 9(b)). The second period ( ) is observed in time interval [ , ] in which the displacement is periodic with period (Fig. 9(c)). The motor rotates at its constant nominal frequency and the meshing stiffness is periodic with a constant period. The dynamic response shows periodic behavior with steady vibration level.

The displacement is non-periodic during the transient regime, who explains that the corresponding spectrum is rich frequency content. The frequency spectrum of the linear displacement of the pinion shaft is presented in Fig. 10. The spectrum indicate the presence of the meshing frequency = 480 Hz and its harmonics × encircled by side-bands. The system eigenfrequencies and are also shown. The presence of side-bands is caused by the transient regime behavior. The high amplitude values of displacement in the transient response are the direct result of the torsional vibration problem caused by the fact that torsional and lateral movements are coupled via the gear mesh [START_REF] Khabou | Study of a spur gear dynamic behavior in transient regime[END_REF].

Conclusions

In this study, the dynamic response of a machine tool system is analyzed in the transient regimes. The variations of input rotational speed are studied. As a consequence of the speed variation, the gear mesh stiffness was modeled with varying period to take into account transient regimes. The torque developed by the electric motor is expressed as a function of the angular position of the crankshaft. The gear mesh stiffness was adequately computed and introduced in the equations of motion of the system and then a numerical simulation was conducted. The numerical results of the dynamic responses of the mechanical system come to confirm a significant effect of the transient regime on the dynamic behavior of a gear set.

  . It is composed by two rotors supported by RBs. The driving rotor (Part 1: P1) has a pinion with teeth, mass and moment of inertia . The driven rotor (Part 2: P2) has a wheel with teeth, mass and moment of inertia .

Fig. 1 .

 1 Fig. 1. Modeling of a single stage spurs gear

Fig. 2 .Fig. 3 .

 23 Fig. 2. Modeling of a single stage spurs gear

Fig. 4 .

 4 Fig. 4. Modelling of Rolling Bearing (RB)

Fig. 6 . 7 .

 67 Fig. 6. Evolution of torque vs time Fig. 7. Mechanical characteristics of the driving and driven systems

Fig. 8 .

 8 Fig. 8. Non-periodic evolution of the meshing stiffness in transient regime

Fig. 12 .

 12 Fig. 12. Evolution in time of the torque applied on the transmission spur gear

Fig. 10 .

 10 Fig. 10. Spectrum of pinion RB1 displacementTo investigate the effect of transient regime on the gear teeth, the inter mesh forces and the transmission error are studied.Fig.11presents the transmission error evolution. It's illustrate that error is time varying and periodic with as the period and its fluctuation is depend of loading conditions.

Fig. 11 .

 11 Fig. 11. Error transmission evolution vs time Fig. 12 indicates the normal evolution of the effort on teeth according to time. It shows that cyclic over loads on teeth can be expected with values reaching 33.10 4 N.

Table 1 .

 1 Parameters

	motors specifications Electric motor type: ABB-MT 90L Electric characteristics 4 poles, 50 Hz, 3 phases 415 V Power (kW) 1.5 Nominal speed (rpm) 1440 Nominal torque (N m) 10 Ratio of starting up torque / 2.7 Ratio of breakdown torque / 3.2 Slip 0.315 Constant of the motor (a) 1.711 Constant of the motor (b) 1.316

Table 2 .

 2 Characteristic parameters of the spur gear

	Teeth numbers Mass (Kg) Inertia moment (Kg•m 2 ) Base circle (m) Module (mm) Rotor torsional stiffness (N•rd/m) Pressure angle (°) Teeth width (mm) Contact ratio	Pinion 20 0.6 2.6 10 -4 45 10 -4 Wheel 40 2.5 0.05 0.11 3 = = 10 5 = 20 23 = 1.6

Table 3 .

 3 The determinate eigenfrequencies of the system

	Eigen frequencies (Hz)
	0 730,5 904,3 1278,3 1424,4 1891,5 2925,6 3481

( ) is the external applied torques vector. It can be expressed as:

( ) = ,0,0,0,0,0,0,-.

Fig. 5. Two dof RB model

Cutting model

The cutting model used in this paper is the dressing and roughing operations. Tangential and normal cutting forces are respectively given by: = ℎ 2.10

with and are terms which reflect the specific cutting pressure, is the width of cut, ℎ is the instantaneous cutting thickness.


Tarik Yahou

James R Dawson

Thierry Schuller

Impact of chamber back pressure on the ignition dynamics of hydrogen enriched premixed flames

Keywords: Spark ignition, Flashback mechanism, Combustion dynamics, Hydrogen flames

The impact of chamber back pressure and H2-enrichment on flashback of premixed CH4-H2-Air flames triggered by ignition dynamics is investigated. Different back pressures are imposed by varying the exhaust flow blockage with perforated plates of different porosities. Experiments measuring the ignition transient, from kernel formation to stable flames and flashback, with increasing levels of H2-enrichment covering a wide range of operating conditions are conducted. Time-series of pressure, velocity and heat release rate are examined to characterise the ignition process. It is shown that the final state of the flame is highly modified by the over pressure caused by the volumetric expansion of the hot gases. Three different scenarios are identified: Soft ignition with direct stabilisation, transient flashback and permanent flashback. The magnitude of the over pressure generated at ignition is significantly amplified by small increases in mean back pressure which, in turn, modifies the flowrate through the injector thereby promoting flashback. In addition to its role as a pressure amplifier, increasing the combustor back pressure also affects the time lag between the over pressure and the peak in heat release rate just after ignition. The impact of H2-enrichment on flashback is investigated by keeping both the bulk injection velocity U b and laminar burning velocity S 0 l constant for flames with 0% to 100% H2-content. It is shown that the magnitude of the over pressure approximately scales with U b /S 0 l and is independent of the total thermal power and H2 enrichment. Although the mean flow conditions and over pressures are kept constant, the stabilisation mode directly depends on H2-content. This shows that flashback triggered by the ignition dynamics cannot be fully described by classical kinematic arguments and a new mechanism is presented.

1 1. Introduction
Given the urgent need to substantially cut CO2 emissions in the face of increasing demand for largescale power generation, the development of carbonfree technologies is essential. Hydrogen fueled gas turbines are ideally placed to play a major role in accelerating the energy transition towards zeroemission power [START_REF] Taamallah | Fuel flexibility, stability and emissions in premixed hydrogen-rich gas turbine combustion: Technology, fundamentals, and numerical simulations[END_REF]. However, the higher flammability, flame speed, and diffusivity of hydrogen and hydrogen-enriched mixtures introduces new challenges that must be addressed [START_REF] Sánchez | Recent advances in understanding of flammability characteristics of hydrogen[END_REF].

In this context, ignition dynamics and flashback raise many issues that need to be considered at the design stage of the engine. For instance, ignition generates an over pressure which can promote flashback inside the injector causing the flame to stabilise in unwanted locations leading to a localised increase in the temperature inside the premixer. Ignition and flashback are phenomena usually studied separately, but the over pressure generated at ignition can lead to flashback.

Spark ignition dynamics have been extensively studied and the process can be categorised into three main steps: i) flame kernel creation, ii) propagation and stabilisation and iii) burner-to-burner flame propagation for multi-injector systems. The fundamental processes underpinning i) and ii) are well documented in terms of spark characteristics, fuel composition and flow velocity [START_REF] Ballal | Ignition and flame quenching in flowing gaseous mixtures[END_REF][START_REF] Mastorakos | Forced ignition of turbulent spray flames[END_REF]. Only recently has research been dedicated to iii), where a large set of data on the light-round process in multi-burner combustors is now available [START_REF] Bourgouin | Ignition dynamics of an annular combustor equipped with multiple swirling injectors[END_REF][START_REF] Machover | Experimental investigation on spark ignition of annular premixed combustors[END_REF]. In single burners, high speed flame visualisations have shown that successful ignition is associated with kernal formation followed by advection and/or propagation to regions close to the anchoring point in order to fully stabilise the combustion reaction. The direction of the instantaneous velocity in the vicinity of the spark has also been shown to play an important role for successful flame growth [START_REF] Marchione | Ignition of turbulent swirling n-heptane spray flames using single and multiple sparks[END_REF]. Recent direct numerical simulations of ignition in high-speed flows above a backward step combustor [START_REF] Pouech | Premixed flame ignition in high-speed flows over a backward facing step[END_REF] have shown that the ability of a flame to propagate into different regions of the flow is essentially controlled by the strain rate in shear regions which can quench the flame and prevent stabilisation.

Compared to constant pressure steady flow combustion systems like gas turbines, a number of investigations on the impact of hydrogen enrichment on the ignition dynamics have been carried out in constant volume chamber or internal combustion engines. Correlations between hydrogen enrichment, the combustion time and the pressure impulse generated at ignition have been established [START_REF] Halter | Characterization of the effects of pressure and hydrogen concentration on laminar burning velocities of methane-hydrogen-air mixtures[END_REF]. For example, it has been shown that even at very low levels of hydrogen concentration (2 Vol.%), the combustion time of hydrogen enriched diesel is significantly reduced without any significant changes in the peak cylinder pressure [START_REF] Liew | An experimental investigation of the combustion process of a heavy-duty diesel engine enriched with h2[END_REF]. Conversely, the effect of hydrogen enrichment has been found to be more pronounced for lean combustible mixtures. In [START_REF] Weinrotter | Laser ignition of ultra-lean methane/hydrogen/air mixtures at high temperature and pressure[END_REF] it was demonstrated that replacing methane with 15 Vol.% of hydrogen led to an increase of the pressure peak by an order of magnitude compared to operating with pure methane. However, it is currently unclear what impact H2 enrichment will have on the ignition dynamics of flames in constant pressure systems.

Furthermore, the interaction between ignition and flashback is not well understood as it has yet to be studied over a wide range of operating conditions and geometries. There have been a number of studies on flashback and their mechanisms leading to the following classification [START_REF] Lieuwen | Burner development and operability issues associated with steady flowing syngas fired combustors[END_REF][START_REF] Kalantari | Boundary layer flashback of non-swirling premixed flames: Mechanisms, fundamental research, and recent advances[END_REF]: (1) core flow flashback, (2) combustion instabilities, (3) vortex breakdown and (4) boundary layer flashback [START_REF] Gruber | Direct numerical simulation of premixed flame boundary layer flashback in turbulent channel flow[END_REF][START_REF] Ebi | Experimental investigation of upstream flame propagation during boundary layer flashback of swirl flames[END_REF]. The second and third mechanisms are not considered in this study. Core flow flashback occurs in freely propagating flames due to a kinematic imbalance between the turbulent flame speed and the local flow velocity [START_REF] Benim | Flashback mechanisms in lean premixed gas turbine combustion[END_REF][START_REF] Lewis | Stability and structure of burner flames[END_REF] and can be avoided by ensuring a sufficient kinematic balance in the burning zone [START_REF] Lieuwen | Burner development and operability issues associated with steady flowing syngas fired combustors[END_REF][START_REF] Benim | Flashback mechanisms in lean premixed gas turbine combustion[END_REF][START_REF] Driscoll | Turbulent premixed combustion: Flamelet structure and its effect on turbulent burning velocities[END_REF]. By contrast, all premixed flames stabilised near walls are potentially susceptible to boundary layer flashback. This can be considered a special case of (1) where the flame propagates into regions of low velocity near the wall but has added complexity such as heat loss. Regardless, these different flashback mechanisms apply to systems operating in quasi-stationary regimes and it is not well understood what impact ignition has on flashback.

In a recent study [START_REF] Prieur | Flame and spray dynamics during the light-round process in an annular system equipped with multiple swirl spray injectors[END_REF], the ignition process was shown to induce a strong modulation in the injector mass flow in both a single sector and annular combustor. This is caused by volumetric expansion of burned gases resulting in a strong pressure perturbation which, in turn, modifies the air flowrate supplied to the burners leading to conditions favourable to flashback. It was observed that during the ignition transient, the flame switched from being anchored inside the injector to a lifted flame. It was later shown that the lower flame branch first flashes back before stabilisation is achieved [START_REF] Töpperwien | Large-eddy simulation of flame dynamics during the ignition of a swirling injector unit and comparison with experiments[END_REF].

The present work aims to improve our understanding of ignition generated flashback for hydrogen enriched flames. It is expected that for increasing hydrogen enrichment, the increased flame speed and diffusivity of hydrogen will increase the susceptibility of the system to flashback compared to premixed gaseous hydrocarbon flames such as methane. To do this, a series of ignition experiments using premixed methane-hydrogen blends with different levels of hydrogen enrichment were performed in a single sector combustor with varying back pressure. First, the impact of the back pressure on the ignition transient is presented and discussed. This is followed by a detailed analysis of two selected operating points to illustrate the effect of different levels of hydrogen enrichment. Finally, the effects of the fuel composition on the ignition dynamics are examined when keeping both the laminar burning velocity and bulk injector velocity constant for different levels of H2-content.

Experimental setup and measurements

The experimental setup is presented in Fig 1a . Different blends of hydrogen, methane and air were premixed upstream and fed into the base of a cylindrical plenum which has a diameter Dp = 65 mm and a length lp = 90 mm. The flow then passes through a converging section before entering a tube with an inner diameter D = 30 mm. Along the centre of the tube is a 6 mm diameter rod which is attached at the top of the plenum and holds a conical bluff body to stabilise the flame that has a diameter D b = 16 mm as shown in Fig. 1b.

The flame is confined inside a cylindrical quartz combustion chamber of inner diameter Dcc = 60 mm. To impose different back pressures on the combustor, a set of perforated plates with varying porosity were used. Each plate had 5 mm holes uniformly distributed in a square pattern as shown in Fig. 1c. The number of holes was varied from N × N = 25 to 16 to modify the porosity σ and vary the pressure drop ∆P at the combustor outlet as shown in Tab. 1. It is important to note that these small changes in the mean static pressure ∆P lead to significant changes in the over pressures generated during ignition.

The inlet velocity was measured 70 mm upstream of the bluff body using a hot wire probe (HW in Fig. 1a). The static and dynamic pressures were measured in the combustion chamber using a Kulite XCS-093-0.35D transducer placed in water-cooled wave 1a. Flame images and time-series of OH* during ignition were measured using a Phantom V2012 highspeed camera with an intensifier (LaVision IRO) and a Hamamatsu H11902-13 photomultiplier (PMT) equipped with a UV filter centred at 310 ± 10 nm. Sampling rates of the PMT and high-speed camera were 52.5 and 10 kHz respectively. It is well established that OH* is a suitable proxy for the fluctuations in the heat release rate in premixed flames.

Care was taken to maintain constant thermal conditions during the ignition experiments. For each test the burner was preheated to a temperature T b = 470 K measured at the center of the bluff-body using a pyrometer. The chamber was then filled with the combustible mixture for 5 s to ensure ignition occurred in a homogeneous mixture. A flame was then ignited using a spark plug positioned on the chamber back plate, approximately 20 mm away from the burner axis, located in the outer recirculation zone of the flow. The ignition system delivered spark energies of 36 mJ at f = 50 Hz. To avoid the saturation of the camera intensifier, a neutral density filter was used to mask the spark from the camera field of view.

The ignition dynamics were investigated for premixed H2/CH4-air mixtures with varying levels of H2 enrichment characterized by the fraction of thermal power defined as PH = PH 2 /(PH 2 + PCH 4 ). Ignition experiments were performed over a wide range of operating conditions from pure methane (PH0) to pure hydrogen (PH100) at different equivalence ratios ϕ for different pressure drops ∆P produced by perforated plates of different porosity placed at the top of the combustion chamber. The experiments were also performed for two bulk flow velocities U b = 2 and 5 m/s at the injector exit under ambient temperature.

Stability maps

Stability maps showing the ignition limits for two bulk flow velocities U b and three different back pressures are shown in Fig. 2. A minimum of five tests were done for each point. Three different ignition scenarios were observed:

• Scenario A: Soft ignition shown in green. Ignition leads to the formation of a flame kernel which propagates through the mixture and stabilises on the bluff-body. Videos of these different sequences are provided as supplemental material. The bottom left region in Fig. 2 withsymbols could not be explored due to the limits of the mass flow controllers. For CH4-air flames, the point at ϕ = 0.6 and U b = 5 m.s -1 also corresponds to the lean blow off limit. The top right region with X symbols is limited by flashback which occurs at lower values of PH as ϕ increases. Note that when ϕ is constant, the thermal power P, remains constant to within ±5% over all cases.

Ub = 2 m/s Ub = 5 m/s X X X X X X X X X X X X X X X X X X - - - - - - - - - - - - - - - 0.
0.8 Plate 16 - X X X X X X X X X X X X X X X X X - - - - - - - - - - - - - - c)
The green shaded region in Fig. 2a shows that atmospheric conditions promote a soft ignition process, Scenario A, for both U b over all ϕ and PH. As there is no flow restriction at the exit of the combustion chamber, the combustion products can freely expand noting that the unsteady back pressure remains in this cases lower than pm ≤ 0.5 kPa during the ignition sequence as will be shown below.

The impact of increasing the back pressure is shown in Figs. 2b and2c. Table 1 indicates that small increases in ∆P from 30 and 140 Pa at the combustor exit measured under stationary cold flow conditions has a strong effect on the ignition process. The Despite the small values of ∆p produced under steady cold flow conditions, they result in significant over pressures pm during ignition causing strong flow modulations through the injector. The amplitude can vary from pm ≃ 2 kPa for PH0 up to pm ≃ 5 kPa for PH100 when ∆P = 100 Pa and from 4 kPa to 8 kPa for PH100 when ∆P = 140 Pa. This modifies the flow conditions promoting the flashback scenarios B and C. These effects are augmented by the higher reactivity of hydrogen leading to a higher acceleration of flame propagation after ignition which, in turn, generates a higher over pressure inside the chamber when the combustor back pressure is increased. This is further investigated in the next section by considering the transient response of the flame and the dynamic pressure.

Ignition dynamics

In the following, two specific flow conditions are considered: a CH4-air flame (PH0) and CH4/H2air flame (PH30) for the same operating conditions U b = 2 m.s -1 and ϕ = 0.7. Figure 2 shows that the PH30 flame switches from Scenario A to B and finally C when the mean back pressure is increased. The PH0 flame is taken as a reference to underline the impact of hydrogen.

PH0

Before discussing the ignition process and transient response, the effect of mean back pressure ∆P on the mean flame shape is considered. Figure 3 shows that the change in back pressure has a minimal effect as all flames show the expected V-flame structure. They feature similar lengths, flame angles and lift-off distances above the bluff body. Slight changes at the end of the flame brush are observed with an increase in ∆P . The different impedance of the plates changes the acoustic reflection which modifies the combustion dynamics at the flame tip. This results in a small change in the mean heat release rate distributions near the top of the flame. Time-resolved pressure, injector velocity, and OH * chemiluminescence were recorded simultaneously during the ignition process. Figure 4 plots their timeseries averaged over five runs and corresponding signal "min-max" envelopes. The instant, t = 0 ms, marks the formation of the first flame kernel which rapidly grows as the fresh gases are quickly consumed leading to a steep increase in OH * reaching peak values between t = 15 -20 ms. It quickly decays to the nominal value I∞ which corresponds to the steady state of the flame. A similar trend in the dynamic pressure is observed but its time lag and amplitude relative to the OH * signal varies with ∆P . To ease interpretation of Fig. 4, the time-series data within the red and gray shaded regions correspond to flashback and ignition transient respectively. The transient ignition period remains roughly constant τp ∼ 35 ms in all three cases.

Without a perforated plate, the ignition process only generates a small over pressure pm ≃ 0.1 kPa, which is not sufficient to modify the flow velocity inside the injector as shown by the black curve in Fig. 4. The addition of Plate 25 and Plate 16 increases the amplitude of pm (blue curves) at ignition and generates correspondingly larger velocity oscillations inside the injector. For Plate 25 the over pressure reaches pm ≃ 2 kPa reducing the flow velocity through the injector to 0.6 m.s -1 at t ≈ 10 ms. As pm decays, the velocity continues to oscillate before mean flow conditions are restored around t = 30 s.

For Plate 16 the over pressure reaches pm ≃ 4 kPa in the bottom plot in Fig. 4 resulting in significant flow oscillations. As the pressure rises, the flowrate inside the injector first decreases but then rapidly accelerates as shown by the large overshoot which peaks around = 5.5 m.s -1 , almost three times above U b , before rapidly decelerating to u ≈ 0.1 m.s -1 . This causes flashback to occur between t = 30 and t = 70 ms highlighted by the red shaded region before stabilising on the bluff body when t > 70 ms. This is shown by following the variation of the OH * signal over the same time period. Although the hot wire cannot detect reverse flow, it is 70 mm upstream, therefore as the velocity approaches zero indicates that reverse flow is likely to occur downstream of the measurement point. Another important feature is the time lag between the over pressure, velocity and OH * timeseries between Plate 25 and Plate 16 which is discussed later in the paper.

PH30

A flame with the same U b and ϕ but a higher hydrogen content, PH30, is now considered. Steady state flame images are shown in Fig. 5 for increasing ∆P . The electrode has been depicted in these pictures for context. It shows that the stability of the flame is very sensitive to the imposed back pressure. Without a plate, an M-shaped flame is stabilized on the bluffbody noting that the flame branch on the left side cannot be seen due to the use of a neutral density filter to avoid saturation of the intensifier during ignition. For Plate 25, flashback occurs and the flame is anchored upstream the bluff-body. In this case the flame exhibits large synchronized self-sustained oscillations leading to a higher OH * intensity in the wake region. Plate 16 leads to permanent flashback as the flame remains stabilised inside the injector. It is worth noting that this behaviour was not observed when flames are ignited with only methane and the H2-content is then slowly increased in a quasi-static way to reach the same enrichment.

The corresponding pressure and OH * times series are plotted in Fig. 6 supported by PIV data in the supplementary material. The hot wire probe could not be used due to flashback. Without any plate, the ignition dynamics leading to soft ignition scenario A is similar to that observed in Fig. 4 for the PH0 case.

The main differences are that the ignition transient (grey shaded region) is slightly shorter τp ∼ 25 ms and pm is approximately twice the value compared with the PH0 case. For Plate 25, the over pressure peak reaches pm ∼ 4 kPa which is of the same order of magnitude as for PH0 with Plate 16 in Fig. 4. In both cases, this leads to Scenario B due to the decrease in the injector velocity is proportional to 1/pm (See the supplementary materials for more details). For Plate 16, the over pressure further increases to pm ∼ 5.5 kPa resulting in permanent flashback of the flame, i.e. Scenario C. It also leads to large pressure oscillations beyond t > 20 ms due to large oscillations of the flame as it cannot properly anchor in the plenum. These changes in the ignition dynamics with H2-enrichment are the consequence of the higher laminar burning velocity which is roughly increased by a factor of 2 from S 0 l = 0.19 m.s -1 for PH0 to S 0 l = 0.35 m.s -1 for PH30 resulting in larger flow accelerations at the exit of the injector as the flame kernel develops [START_REF] Prieur | Flame and spray dynamics during the light-round process in an annular system equipped with multiple swirl spray injectors[END_REF].

The results for PH0 and PH30 show that the main mechanism driving flashback is the transient, high amplitude over pressures. The over pressure amplitude is very sensitive to the mean combustor back pressure as shown by the addition of the porous plates. During the ignition transient, the flame acts as a piston moving at the absolute flame speed Sa = w • n [START_REF] Poinsot | Theoretical and numerical combustion[END_REF] which pushes away fresh gases [START_REF] Töpperwien | Large-eddy simulation of flame dynamics during the ignition of a swirling injector unit and comparison with experiments[END_REF]. The addition of porous plates lead to small increases of the mean back pressure amplifying the amplitude of pressure perturbations due to ignition, which in turn modify the flow response through the injector promoting flashback. However, the results also suggest that the synchronization of the pressure and OH * signals also plays an important role in triggering flashback. This makes it difficult to isolate the role of over pressure and H2content on flashback.

Impact of hydrogen enrichment

An attempt is made to isolate effects of combustion properties of hydrogen from the flow dynamics by keeping both the bulk velocity and the laminar burning velocity fixed at U b = 5 m.s -1 and S 0 l = 0.25 m.s -1 respectively for increasing levels of hydrogen enrichment. To maintain a fixed burning velocity S 0 l , ϕ was varied in accordance with Tab. 2. Keeping U b and S 0 l constant necessarily results in a change in the thermal power.

Figure 7 shows average time-series of dynamic pressure measured in the combustor and OH * emission during ignition for different levels of H2enrichment (PH%). Averages were obtained for a minimum of 5 runs for each PH. For brevity, the results are only presented for case Plate 25. Figure 7a shows that keeping U b and S 0 l constant scales the response of the time-series. All levels of PH achieve approximately the same over pressure pm ≃ 5 kPa, which consequently leads to a similar flow perturbation within the injector. This result confirms that for a constant laminar burning velocity S 0 l , the effects of thermal expansion and flame propagation during ignition remain comparable and are independent of the level of H2-enrichment.

Although good collapse of the over pressures with S 0 l is found, Fig. 7b shows that the flame stabilization process differs with PH. From PH0-PH20 soft ignition and stabilization processes occur with Scenario A. At PH40, temporary flashback occurs with Scenario B between 15 ms ≤ t ≤ 40 ms before stabilizing at the bluff body. Once PH 2 > 40% all the flames permanently flashback with Scenario C immediately after ignition. This shows that flashback cannot be a function of S 0 l only and that the level of H2enrichment needs also to be considered.

The effect of H2-enrichment on the flame response cannot be fully explained using classical kinematic arguments which suggest that flashback is mainly driven by an imbalance in the kinematic equilibrium between the local flow and burning velocities. For example, simple flashback criteria are often defined when the ratio (U b /S 0 l ) drops below a critical threshold value (U b /S 0 l )c [START_REF] Kalantari | Boundary layer flashback of non-swirling premixed flames: Mechanisms, fundamental research, and recent advances[END_REF][START_REF] Benim | Flashback mechanisms in lean premixed gas turbine combustion[END_REF]. In these experiments, the ratio (U b /S 0 l ) is kept constant, yet the susceptibility to flashback increases with increasing hydrogen enrichment. Furthermore, since the over pressure remains also the same in all cases, the increased propensity for flashback cannot be due to differences in the injector flow dynamics. Preferential diffusion effects must therefore play a crucial role. One also needs to keep in mind that the unsteady flame front is highly wrinkled and propagates anisotropically. These local effects cannot be captured using the photomultiplier (PMT in Fig. 1).

Two conditions need to be simultaneously met to allow the flame to propagate into the injector. First, the ratio (U b /S 0 l ) needs to drop below the critical flashback limit (U b /S 0 l )c, second, the flame must be within the vicinity of the injector. Synchronisation of these two events is therefore essential. The first condition is set by the amplitude of the over pressure and the critical threshold (U b /S 0 l )c which is kept constant. The critical threshold is expected to be lower for H2-enriched flames given that S 0 l is lower than the real flame displacement speed S d . For lean hydrogen mixtures, S d exceeds S 0 l due to the imbalance between heat and species diffusion in flames characterized by non-unity Lewis numbers which generates reaction layer instabilities [START_REF] Berger | Intrinsic instabilities in premixed hydrogen flames: parametric variation of pressure, equivalence ratio, and temperature. part 2-non-linear regime and flame speed enhancement[END_REF].

In the cases explored during this study, pressure peaks larger than pm ≃ 4 kPa led to a temporary reverse flow inside the annular injector as highlighted in Fig. 8.a. This figure plots the non-filtered time series of a single run for PH0 flame in Tab. 2. It clearly shows that the pressure perturbation leads to a temporary reverse flow between t = 8 ms and t = 13 as the velocity signal measured by a hot wire (HW in Fig. 1) first drops to zero at t = 8 ms, then rapidly increases before dropping down again to zero at t = 13 ms followed by a large overshoot reaching a peak around u ≃ 14 m.s -1 . Since the hot wire measures only the the absolute component of the velocity, the real velocity signal is reconstructed and plotted in Fig. 8.b with negative velocities between 8 ms ≤ t ≤ 13 ms that drops down to u ≃ -6 m.s -1 at t ≃ 11.6 ms. This instant approximately matches the instant t = 11.2 ms when the pressure reaches its maximum. This observation indicates that the first condition for flashback is always satisfied in Fig. 7 regardless of the H2-content in the fuel. Despite the strong reverse flow, the flame smoothly ignites following ignition Scenario A in Fig. 8 demonstrating that the first condition based on the kinematic balance between the reactant flow velocity and burning velocity is a necessary but not sufficient condition to trigger flashback.

The second condition considers the location of flame propagation and advection in the volume and cannot be verified via PMT data. Nevertheless, the synchronization of these two events can be investigated by evaluating the time-lag τ between the peak over pressure and OH * time-series as shown in Fig. 9. It was previously shown that the instant corresponding to the peak pm also coincides with the minimum velocity inside the injector (see Fig. 8). If the OH * peak is considered to mark the instant when all the fresh gases have been consumed and therefore, corresponds to the time when the flame is stabilized on the bluff-body, then, as τ tends to zero, the time synchronisation condition is satisfied increasing the risk of flashback. In particular, if the velocity inside the injector is sufficiently low, both conditions are met and the lower flame branch can propagate upstream into the injector causing flashback. On the other hand, when τ is large enough, the sudden velocity drop inside the injector occurs before the flame has time to reach the injector orifice and flashback is avoided.

Based on the data, a critical value τc ∼ 1 ms is identified as the minimum time lag τ > τc above which the synchronisation condition is not satisfied and the ignition leads to Scenario A (green symbols in Fig. 9). Conversely, when τ < τc, the flame flashes back and leads to Scenario C shown by the red symbols. When τ ∼ τc, the ignition develops Scenario B (orange markers) characterized by transient flashback. Figure 9 also shows that the delay τ decreases with increasing ∆P . This decay is partly related to the inertial time τ f needed to displace the column of fluid by thermal expansion of the burned gases which can be considered as the sum of two time scales τ f = τ l + τ δ . The first delay τ l = lc/Sa denotes the time required for the flame propagating at absolute speed Sa to cross the length lc of the combustion chamber. The additional delay τ δ = δσ/Sa is due to the inertia of the fluid exhausting through the holes of the perforated plate, where δ is the sum of the plate thickness and an end correction which depends only on the diameter of the holes and σ is the plate porosity. In the conditions explored, τ l is fixed and since τ δ is proportional to σ, τ f decreases with increasing ∆P .

For a fixed porosity, τ also reduces in Fig. 9 with hydrogen enrichment. This might be linked to the drop of the ignition time delay τi, defined here as the time interval between the first spark and the sudden rise of OH * emission. With this definition, experiments show that τi decreases from τi ∼ 5 ms for PH0 to τi ∼ 0.2 ms for PH100. However, this gap does not correspond to those observed in Fig. 9 when one switches from PH0 to PH100 for a given back pressure. A drop in τi also implies an increase in the flame speed. This would lead to an increase in over pressure as shown in Section 4, which is not observed here.

Resistance to stretch can also affect flame propagation during ignition. Pouech et al. [START_REF] Pouech | Premixed flame ignition in high-speed flows over a backward facing step[END_REF] showed that the development of a flame kernel initiated in a recircu- lation zone can be quenched when entering in a highshear region of the flow where the local strain rate exceeds the critical extinction strain rate. In these experiments, the shear and strain rates produced by the annular jet in cold flow conditions are the same for all combustible mixtures. Due to their higher resistance to strain, flames with higher levels of H2-enrichment can propagate into regions of more intense shear and propagate through the annular jet shear layer.

This interpretation is qualitatively supported by Fig. 10 which shows instantaneous OH * images between t = 3 ms and t = 6 ms after ignition. For the PH0 flame at the top in Fig. 10, the heat release rate distribution is not uniform along the flame front. The upper flame branch is distinguished by a high intensity and the lower flame branch, slightly lifted, with a lower light intensity. The overhead views show that the flame cannot cross the jet shear layer and stabilises along the shear layer by propagating around it. The PH100 flame at the bottom in Fig. 10 shows a homogeneous heat release rate distribution along the flame front with a lower flame branch that is already propagating in the upstream direction through the injector for t ≥ 5 ms. The overhead view shows the flame propagating through the jet shear layer towards the bluff-body which results in flashback when the flow rate drops. Finally, for lean hydrogen enriched flames thermo-diffusive effects can counteract the effects of strain enabling flame propagation through regions of high shear due the effects of preferential diffusion which leads to local enrichment along the reaction front and a higher local burning velocity [START_REF] Berger | Intrinsic instabilities in premixed hydrogen flames: parametric variation of pressure, equivalence ratio, and temperature. part 2-non-linear regime and flame speed enhancement[END_REF].

Conclusion

The impact of chamber back pressure and increasing H2-enrichment on the ignition dynamics of premixed bluff body flames has been investigated. It has been shown that the propensity of the flame to flashback depends on both the amplitude of the over pressure produced by ignition and the level of H2-enrichment. Using perforated plates of different porosity to vary the back pressure, it was shown that small changes in the pressure drop amplify the pressure perturbations at ignition leading to very large over pressure amplitudes and promoting flashback. The large amplitude over pressures reduce the reactant flowrate at the injector exit facilitating flame propagation in the upstream direction.

The amplitude of the over pressure was shown to scale with the ratio of bulk velocity to laminar burning velocity. However, for a fixed over pressure the risk of flashback was found to increase with increasing levels of H2-enrichment. Two mechanisms were considered to explain this behavior, including a shift of the quasi-static flashback limit due to an increase of the flame speed and an increased resistance to stretch with a higher extinction limit of the propagating reaction layer. The latter effect is due to the high diffusion rate of hydrogen which enables flames to propagate through regions of high shear whereas mixtures with little or no hydrogen are extinguished.

For a given over pressure, the risk of flashback increases when the peak of the over pressure coincides with the peak of heat release rate. This synchronization is necessary to get a flame front close to the burner outlet when the flow velocity drops inside the burner. When the time lag between the peak pressure and heat release rate increases, the risk of flashback is substantially reduced.

Fig. 1 :

 1 Fig. 1: (a) Experimental setup and relevant dimensions in millimeters. (b) Zoom of the burner outlet. (c) Combustion chamber exhaust with plate 25: N × N = 25 holes.

Fig. 2 :

 2 Fig. 2: Stability maps for (a) W/o plate , (b) Plate 25 and (c) Plate 16 configuration. The green shaded region represents the ignition Scenario A, the orange one Scenario B and the red one Scenario C.

Fig. 3 :

 3 Fig. 3: Mean flame images for PH0 at ϕ = 0.7 and U b = 2 m.s -1 for three back pressures presented in Tab. 1.

Fig. 4 :

 4 Fig. 4: Time traces of chamber pressure, normalized OH * light intensity and local flow velocity for PH0 at ϕ = 0.7 and U b = 2 m.s -1 . The time traces are averaged over 5 runs. The range of pressures varies in each plot.

Fig. 5 :

 5 Fig. 5: Mean flame images for PH30 flame at ϕ = 0.7 and U b = 2 m.s -1 .

16 Fig. 6 :

 166 Fig. 6: Time traces of the chamber pressure and normalized light intensity averaged for PH30 flame at ϕ = 0.7 and U b = 2 m.s -1 . Solid lines stand for the average signal over the runs. The envelope of all runs is also plotted.

Fig. 7 :

 7 Fig. 7: Chamber pressure (top) and light intensity (bottom) signals for the six flames presented in Tab. 2 with Plate 25.

Fig. 8 :

 8 Fig. 8: Time traces of the chamber pressure (blue line), OH * light intensity (red line) and injector flow velocity (black line) for a single run for PH0 flame in Tab. 2. a) shows the absolute injector flow velocity measured by the hot wire and b) denotes the real velocity signal.

Fig. 9 :

 9 Fig. 9: Time lag between the chamber pressure and light intensity peaks for the three configurations presented in Tab. 1. Green symbols: ignition Scenario A. Orange symbols: Scenario B. Red symbols: Scenario C. The black dashed line corresponds to τc = 1 ms.

Fig. 10 :

 10 Fig. 10: Instantaneous OH * snapshots highlighting flame propagation after ignition for S 0 l = 0.25 m.s -1 and U b = 5 m.s -1 . Top row: PH0. Bottom: PH100. Operating conditions are given in Tab. 2. On the top right of each image, a sketch of the overhead view is drawn to delineate the flame position represented in blue.

Table 1 :

 1 Plate characteristics, plate porosity σ and pressure drop ∆P during cold-flow stationary conditions. The bulk velocity U b is defined at the location of the hot wire probe.

		N × N	σ	∆P 2 m.s -1 5 m.s -1
	W/o Plate	-	1	∼ 0 Pa 30 Pa
	Plate 25	5 × 5	0.17	25 Pa 100 Pa
	Plate 16	4 × 4	0.11	40 Pa 140 Pa
	guide, labelled M1 in Fig.

Table 2 :

 2 Operating conditions for different levels of H 2enrichment keeping both the bulk velocity U b = 5 m.s -1 and the laminar burning velocity S 0 l = 0.25 m.s -1 fixed.

	Cases	ϕ	PCH 4	PH 2	P [kW]
	PH0	0.78 100%	0%	8.4
	PH20	0.68	80%	20%	7.6
	PH40	0.60	60%	40%	6.9
	PH60	0.53	40%	60%	6.3
	PH80	0.48	20%	80%	6.0
	PH100 0.43	0%	100%	5.7

Acknowledgments

The authors acknowledge support from the NCCS Centre, funded under the Norwegian research program, Centres for Environment-friendly Energy Research (FME) (Grant 257579/E20).

Supplementary material

Videos of the three different scenarios, effect of the chamber length on the ignition dynamic as well as instantaneous PIV snapshots for PH30 flame at U b = 2 m.s -1 and ϕ = 0.7 are provided as supplementary material.


Marc Berthillier

Christian Dupont

Ritwik Mondal

Jean-Jacques J Barrau

M Berth Ii Lier

R Mondal Blades

Blades Forced Response Analysis With Friction Dampers

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Friction dampers are widely used in axial turbomachinery components to reduce vibration amplitudes and therefore to increase high frequency fatigue life of blades and disks. The analysis of such systems is difficult due to the nonlinear nature of the friction law. Nevertheless, fr om an industrial point of view, it becomes necessary to be able to conduct analysis in order to optimize friction dampers.

Numerical integration methods, although very effective to handle multiple sticking and complicated damping laws, lead to prohibitive computational effort if parametric studies and frequency response curves are to be obtained within the design process. If one is only interested in the stationary behaviour under periodic loading, frequency domain methods are more suited.

The harmonic balance method (Nayfeh and Mook [START_REF] Nayfeh | Nonlinear Oscillations[END_REF]) has been applied to study friction damping in turbine blades by Griffin [START_REF] Griffi N | Friction Damping of Resonant Stresses in Gas Tur bine Engine Airfoils[END_REF]. Although computationally efficient, this method assumes the response to be limited to a single harmonic approxi mation. The results are adequate for continuous slip motion, but not for stick-slip motion where higher harmonics become important.

Multiharmonics methods have been introduced by various researchers. Lau and Cheung [START_REF] Lau | Amplitude Incremental Variational Principal for Nonlinear Vibration of Elastic System[END_REF]. Lau et al. [START_REF] Lau | Variable Parameter Incrementation Metbod for Dynamic Instability of Linear and Non linear System[END_REF][START_REF] Lau | Incremental Harmonic Balance Method With Multiple Time Scales for Nonlinear Aperiodic Vibrations[END_REF], proposed an incremental harmonic balance (IHB) method to solve nonlinear periodic and quasi periodic vibrations of continuous systems.

Pierre, Ferri and DoweU [START_REF] Pierre | Multi-Harmonic Analy sis of Dry Friction Damped Systems Using an Incremental Harmonic Balance Method[END_REF] extended the JHB method to ana lyse the steady state response of systems with dry friction. Ferri [START_REF] Ferri | On the Equivalence of tbe Incremental Harmonic Balance Method and the Harmonic Balance-Newton Raphson Method[END_REF] showed the equivalence between the IHB method and the Galerkin/Newton-Raphson (GNR) method proposed by Urabe [START_REF] Urabe | Galerkin' s Procedure for Nonlinear Periodic Systems[END_REF] to analyse non linear periodic systems. More recently, sev eral authors as Ling and Wu (9), Cameron and Griffin [JO]. Cardona et al. [START_REF] Cardona | A Multihar monic Method for Non-Linear Vibration Analysis[END_REF]. improved the numerical efficiency of the IHB/GNR method by the introduction of the FFT algorithm. However, concerning the application of these methods to industrial cases, with many degrees of freedom, very few work has been published.

The purpose of this paper is to apply the GNR method with FFT algorithm to predict the forced response of turbomachinery blades with friction dampers in an efficient way. The size reduc tion of the system to be solved has been performed in two steps, first by a Craig-Bampton component mode synthesis and second by a condensation technique in the frequency domain on the non linear degrees of freedom only. The accuracy of the method is discussed and the results are compared with experiments conducted on a cantilever beam with a dry friction damper attached. Finally, the application of the method to a model fan blade is presented.

Theory

Finite element models of blades usually comprise a great number of degrees of freedom ( dof). thus to simulate forced response with friction dampers, one has first to tackle the prob lem of model reduction. Dowell [ 12) proposed a component mode analysis with the modes of the structure without dampers. In order to handle systems of smaller size, we prefer to use a Craig-Bampton component mode synthesis, where the retained dof include those where frictional forces are applied. In this paper, those dof are called nonlinear dof. The set of retained dof may be restricted to the only non linear ones, but may include others where one wishes to know the displacement with out using a restitution procedure. The influence on the results of the r etained dof choice will be discussed later. The condensed mass and stiffness matrices are given directly by the finite ele ment code after component mode synthesis.

The forced response with friction is computed with a GNR method using FFT algorithm. The system of equations of motion can be represented by:

l(x(t)) -f(t) -g(x, x) = 0 ( I)
where L( . ) is the differential operator representing the linear term, f(t) is the external forcing te1m, g(x, x) represents the nonlinear friction force and x(I) is the system response. All are vector functions.

It has been shown by Korkmaz et al [START_REF] Korkmaz | Theoreti cal Dynamic Analysis of a Cantilever Beam Damped by a Dry Friction Damper[END_REF], among others, that for a friction damper, stick-slip mechanism is the most efficient in reducing amplitude. As a classical friction Jaw is able to reproduce stick-slip well, it has been used throughout this work.

Thus when sliding occurs between the two dof i and j. the fri ction force acting on the dof i is:

g'(x, xY) = -8d sgn i:if if .xu * 0
where xii = t1ij, and 8d is given by a coefficient of friction times normal force to the sliding surfaces.

For numeri cal purposes, the friction law g has been approxi mated by: = -8d sgn j;V if I xlil > ll.

where A is a numerical parameter.

If we assume the solution to be:

x(t) = � (! X co + i Xc,, cos nwt + Xs11 sin nwt ) (2) 7r 2 11�1
the Galerkin procedure leads to the spectral form of (1 ) (see references 9, 10 and 11, for more details) :

A• X -F -G(X) = 0 (3)
where the matrix A and the vector F are easily obtained from L and f. X is the vector given by:

The vector G(X) is given by: where and GT(X) = (G�0, G�i. G'J' i. ... , G��. G'{N) r2ir lw

Gc .. (X) = Jo cos (nwt)g(x)dt f2ir lw

Gs,,(X) = 0 sin (nwt)g(x)dt

The solution vector X of Eq. ( 3) is found by an iterative Newton-Raphson procedure. If we write the residual vector R(X) as:

R(X) = A•X-F-G(X)
the iterative system to be solved is:

J(X)•!:.X = -R(X) (4)
where !:.X is the increment of the displacement vector X to be found and J ( X) is the jacobian matrix:

J(X) = 8R =A -

ao ax ax

In order to calculate G(X) and 8G/8X by numerical FFT, the functions are discretised with NP points per period (NP > 2N, N = number of harmonics) :

!:.t = � wNP xk = x(k!:.t) k from 0 to NP-1 ik = i(k!:.t) j ( ,ij ) 1 ,ij ' f 8 k xk = -gd-xk 1 A. I .x f I > A
Each component i of Gc11(X) and Gs11(X) is approximated by:

G�n(X) = L cos _!!__ nk gi(x f ) NP-1 ( 2 )
k�o NP These terms are zero if the dof i is not affected by friction.

The matrix 8G/ 8X is composed of (2N + 1) x (2N + 1) block matrices of dimensions p x p where p is the number of dof of the initial temporal system. 8Gc11I 8Xcm is the p x p block corresponding to the line c n and to the column cm. The other generic blocks are: BGs,,I 8X,m, 8G,n18Xcm• 8Gcn18X,m.

The only non-zero terms of these matrices correspond to those affected by friction. For example, the generic term i, j of the first block matrix is numerically calculated by: en

I cos -nk --.- 8Gj ex) NP -L (27r ) 8g' k •c.x' k u) 8x fm k�o NP 8X {,m
The expression of the differential terms is not straightforward, it is obtained by: 8gt = 8gt. a.x f 8x {,,, a.x f 8X !,,,

Since:

dxk = -L -n srn -nk d(X e n -X {,,,) ,ij W 2 ( N , ( 27r ) j .
7r n�1

NP we can write:

( 27r ) . . )

+ n cos NP nk d(X � . . -X f11) 8x f mw2 ( • ( 27r k )) --. =--sm -m . 8X {,m 7r
NP So, we obtain:

8G�,, _ mw 2 ( N �i ( 27r k ) . ( 27r . k ) 8gi) --. --- "°' cos -n sm -m -.. 8X {,m 7r M NP NP 8.X f = mw _.::'.._ I, sin _!!__ k(n + m) � ( NP-I ( 2 ) 8 j 27r M NP 8 .Xf --L sm -k(n -m) -. . W NP-1 , ( 27r ) 8gi) 27r M NP 8 .Xf
In prac�ice, we begin by numerically calculating the terms 8gi/8x% for k varying from 0 to NP-1, with the formulas:

1 8g� = 8 .x f -g d � =0 if if l x f I ::5 A. l x f I > A.
The terms 8G�11/ 8X 1m can then be calculated, for n and m varying from 1 to N, by using an FFT algorithm.

In order to reduce even more the computational time, the iterative system ( 4) is reduced to the non linear dof only, i.e. the dof affected by friction. The solution vector X can be partitioned according to the linear and non linear dof as X = ( X 1 11 , X111) , and the system ( 4) can be written:

It follows immediately that it is equivalent to solve one system of smaller size: and to restore the linear dof by: !:.X1 n = J b, 1 (R1 n -11 11,111!:.X .. 1 )

This reduction technique is specially interesting when the number of normal modes and retained linear dof of the compo nent mode synthesis is important.

Experiment

In order to gain a better understanding of friction damping and to have experimental data to validate numerical methods, a simple experiment has been realized for a cantilever beam excited by an electrodynamic shaker and damped by a dry fric tion damper. A brief description of the experimental set up and of the results will be given here. The reader is referred to Kork maz et al. [START_REF] Korkmaz | Analyse Dynamique Experimentale d'une Poutre Amortie par un Amortisseur a Frottement Sec[END_REF] for more details. The design of the friction dampers is shown in Fig. 2. A friction damper is applied against each of the two plates fixed to the beam, in order to preserve symmetry. It is possible to apply a variable normal load on the contact area of each damper.

An electrodynamic shaker is attached to the beam at a dis tance of 1.28 m from the clamped end.

Accelerometers are placed on the beam, perpendicular to span, at the friction damper (point 2) and at the shaker (point 1) locations. Measurements. A normal load FN is applied on the fric tion dampers. Then the beam is excited with a periodic force P cos (wt), whose frequency is swept over a given range.

Measurements have been done for several levels of excitation force P and normal loads FN. For all measurements, the transfer functions between the excitation force P and the responses at The beam tested previously has been chosen for the validation of the computations.

Validation has been conducted as follows:

-Creation of the Finite Element Model.

-Creation of 6 reduced models (1 or 2 retained dof, 1, 3 or 10 normal modes), which have been compared to the Finite Element Model with respect to eigen-frequencies and mode shapes predictions.

-Analysis of the multiharmonic, friction damped, forced re sponse. The influence of the reduction technique and the number of harmonics has been studied.

-Analysis of the temporal, friction damped, forced response. Component Mode Synthesis. The choice of the retained dof is an important step in the model reduction. The smallest set of retained dof is composed of those where friction dampers apply. In our beam case, this set is reduced to a single translation dof and is called setl. For large structures, to avoid a lengthy displacement restitution procedure, one may wish to retain addi tional dof where the displacement is to be known. For our beam, we add to set l the displacement at the shaker location, to form the set2. We performed several component mode syntheses (CMS), with 1, 3 and 10 normal modes (NM), for setl and set2.

Finite Element

The 6 reduced models obtained, have been compared to the initial Finite Element Model for two eigenvalue analyses. One with a spring of value equal to 2.4 x 10 7 N Im attached at the damper location and the other with the displacement at the same location free to move. These configurations are called hereafter respectively fixed and free boundary conditions. The eigen frequencies obtained with the various models are presented on table 1. The modal displacement ratio (u2/ u3) X 100 for these models is reported on table 2, where u2 and u3 are respectively the displacements perpendicular to span at the damper location and at the free end of the beam.

To analyse the first beam mode forced response, one may think that the 1 NM models are pertinent. The 10 NM setl Harmonic Analysis. We chose to simulate forced response with friction damping for the first mode. This mode has been chosen because we have experimental results. The results ob tained with the reduced models presented previously, will be compared with one another, in order to evaluate the effective ness of the component mode synthesis technique, combined with friction modeling. Then the influence of the number of harmonics in the response will be studied.

We analyse the beam in the configuration that leads to the results presented on Fig. 4 for the excitation force levels P of 0.938 N, 25 N and 38.438 N. When we simulate the forced response with a single harmonic, all models give equivalent results. The values of u 1 IP X 105 over the frequency range of 20 Hz to 40 Hz are plotted on Fig. 5 and the maxima are reported in table 3. On Fig. 5, the results of all models superimposed the reference model (10 NM set 1) results, except for the 1 NM set 1 model. In fact, we saw that this model does not represent perfectly the modal displacement at the damper location for the free and fixed boundary conditions (see Table 2). In the follow- ing part of this work, the 1 NM set 2 model will be used to study the fi rst mode forced response.

We can say that component modes synthesis is an effi cient way to reduce the initial Finite Element model for studying forced response with dry friction. Nevertheless it is necessary to check the reduced models ability to represent the eigen fre quencies and mode shapes of the initial model for the free and fi xed boundary conditions. We find that, for our case studies, small size models are sufficient.

The influence of the harmonics number has been studied for the first mode. The results are presented on Table 4 and on Figs. 6 and7. We see that the number of harmonics has a great influence on the response when stick-slip occurs. We find for P = 25 N a 35 percent difference on the maximum response amplitude when computing the response with a single harmonic or with seven harmonics.

Finally the frequency response functions have been computed for all experimental excitation force levels with seven or nine harmonics. These curves are presented Fig. 8 and are to be compared with those of experimental results shown in Fig. 4. We can say that the results agree quite well. The computational time (CPU time) is much more dependent on the system size to be solved than on the harmonics number. We fi nd, as an average, that the resolution of a system with 3 equations needs 30 percent of the CPU time required for 11 equations, whereas a one harmonic analysis needs 70 percent of the CPU time required for a seven harmonics analysis.

Temporal Analysis. Temporal analyses have been con ducted with the initial Finite Element model. The steady state responses have been compared with the experimental and multi harmonics (GNR) method results. The computed values are presented on Table 5.

We can first note that the GNR method fits very well with Time Integration method except for low excitation forces. How ever these low excitation responses do not correspond to the highest effi ciency of friction dampers and are of minor interest from a design point of view. The correlation with experimental results is acceptable and the discrepancies are due to modelisa tion, not to numerical methods.

Application

The method is now applied to a model fan blade, damped by a blade to ground damper and excited by a concentrated periodic force. The blade root is clamped and the first flexion Results. Figure 10 shows, for a fixed excitation force, the displacement amplitude at the blade tip as a function of excita tion frequency for different normal forces (from 0 to Fl) ap plied at the frictional damper. A single harmonic approximation for the response has been considered in these computations. We note that for the normal forces Fl and 0, the frequencies at which the amplitudes are maximum correspond to the first eigen frequencies of the blade, for fixed and free conditions at damper location. Numerical results correlate qualitatively well with using component modes synthesis and to perform the eigen frequencies analysis as well as the nonlinear temporal computa tions.

Conclusions

A numerical method to study the forced response of blades with friction dampers has been presented. This method com bines a multiharmonic analysis with a Craig-Hampton compo nent mode synthesis and is easy to use for industrial analysis as presented here for a fan blade. It allows damper optimization within the design process.

The accuracy of the method has been tested on a cantilever beam. It has been shown that, with component mode synthesis, small size reduced models are able to represent properly the structure dynamics for friction damping analysis. Higher har monics of the response greatly influence the result of the GNR method, specially when stick-slip occurs. A 35 percent differ ence has been found, in one case, on the displacement amplitude between a one and a seven harmonics analysis. As the highest effectiveness of a friction damper is obtained when stick-slip is present, multiharmonic analysis is of great interest if one wants to optimize dampers. In most of the useful configurations, the analysis with seven harmonics fits reasonably well with temporal analysis and with experimental results.

Fig. 1

 1 Fig. 1 Beam and dampers overview

  Experimental set up. The beam, represented in Fig. l, is a straight iron tube of 1.33 m length, with a rectangular cross section of 0.08 m by 0.04 m and a thickness of 0.002 m. Two very stiff plates are fixed to the beam at a distance of 0.318 m from the clamped end.

Fig. 2 .Fig. 3

 23 Fig.2Friction damper mechanism

point 1 ( 8 NFig. 4

 184 Fig. 4 Experimental frequency response functions for the first beam mode

  Model. A finite element model of the beam with I 0 beam elements has been done. Additional lumped pa rameters have been added to represent the mass ( 2 kg) and inertia ( 8.5 X 10-4 kg.m2) of the added plates and the moving mass of the shaker (0.133 kg). The stiffness of the friction dampers is equal to 2.4 x 107 Nim. The material properties of the beam are: Young's modulus E = 2.0 X 1011 Pa, Poisson's ratio v = 0.3 and mass density p = 7800 kg/m 3• The sliding and sticking eigen-frequencies of the model have been found to be respectively 25.4 Hz and 35.9 Hz. The agreement with the experimental results is acceptable. Structural damping for the first mode was estimated to be e = 1.95 per cent when the damper is sticking (35.9 Hz) and e = 2.55 per cent when the damper is sliding (25.4 Hz). Rayleigh damping was then assumed in the form a K + {3 M where a = 2.6 x 1 o-s and {3 = 7.49.

Fig. 5

 5 Fig. 5 Model influence on the first beam mode forced response

x 1 Fig. 6 o 3 5 Fig. 7 Fig. 8

 163578 Fig. 6 Influence of the harmonics number on the first beam mode forced response P = 38.438 N

Fig. 9

 9 Fig. 9 Fan blade FEM component tests in progress. Detailed comparison of results will be reported elsewhere.Throughout this work, the general purpose Finite Element software code SAMCEF was used to obtain the reduced models

  Fig. 10 Numerical frequency response function for the first blade mode

Table 1

 1 Eigenfrequencies (Hz); Mode 1-damper location free or fixed

			free	fixed
	Full FEM		25.40	35.91
		1NM	25.40	35.91
	CMS set 1	3NM	25.40	35.91
		10NM	25.40	35.91
		1NM	25.40	35.92
	CMS set 2	3NM	25.40	35.92
		10NM	25.40	35.92

Table 2

 2 Displacement ratio u2/u3 x 100; Mode 1-damper location

	free or fixed
			free	fixed
	Full FEM		8.982	1.440
		1NM	8.970	1.439
	CMS set 1	3NM	8.982	1.440
		10NM	8.982	1.440
		1NM	8.981	1.439
	CMS set 2	3NM	8.982	1.440
		10NM	8.982	1.440

Table 3

 3 Model influence on the first beam mode forced response Max(u1/P x 10 5 ) m/N over 20Hz-40Hz

	p		CMS set 1		CMS
					set 2
		1 NM	3NM	10NM	1NM
	0.938N	38.2	38.2	38.2	38.5
	25.0N	4.9	5.1	5.1	5.1
	38.438N	11.4	12.2	12.2	12.1

Table 4

 4 Influence of the harmonics number on the first beam mode

	forced response
	Max(u1/P x 105 ) m/N over 20Hz-40Hz
	p	0.938N	25.0N	38.438N
	1H	38.46	5.07	12.11
	3H	38.46	7.57	13.67
	5H	38.46	7.24	14.83
	7H	38.46	6.84	14.95

Table 5

 5 Analysis and experiment comparison for the first beam mode

	forced response
	5 ) m/N over 20Hz-40Hz GNR Max(u1/P x 10 Time lnteg. Exp. p 0.938N 38.46 36.70 36.10
	25.0N	6.84	6.84	7.89
	34.375N	11.45	11.49	9.31
	38.438N	14.95	14.90	14.00