Swann Marx
email: swann.marx@ens-cachan.fr.

Eduardo Cerpa
email: eduardo.cerpa@usm.cl.

Output Feedback Control of the Linear Korteweg-de Vries Equation*

This paper presents the design of an output feedback control for a linear Korteweg-de Vries equation. This design is based on the backstepping method which uses a Volterra transformation. An appropriate observer is introduced and the exponential stability of the closed-loop system is proven.

I. INTRODUCTION

The Korteweg-de Vries equation, introduced in 1895 by the Dutch mathematicians Diederik J. Korteweg and his student Gustav de Vries, describes approximatively the behavior of long waves in a water channel of relatively shallow depth. This nonlinear partial differential equation, described by u t (t, x) + u x (t, x) + u xxx (t, x) + u(t, x)u x (t, x) = 0, [START_REF] Cerpa | Control of a Korteweg-de Vries equation: a tutorial[END_REF] has been deeply studied in the controllability sense (see [START_REF] Rosier | Control and stabilization of the Kortewegde Vries equation: recent progresses[END_REF], [START_REF] Cerpa | Boundary controllability of the Korteweg-de Vries equation on a bounded domain[END_REF], [START_REF] Cerpa | Control of a Korteweg-de Vries equation: a tutorial[END_REF] and the references therein). By considering different boundary conditions on an interval [0, L] and different boundary actuators, we get control results of different nature. Roughly speaking, the system is exactly controllable when the control acts from the right endpoint x = L, and nullcontrollable when the control acts from the left endpoint x = 0.

In this article we focus, as a first step opening further research, on the stabilizability problem for the linear Kortewegde Vries equation with a control acting on the left Dirichlet boundary condition. The studied system can be written as follows:

     u t + u x + u xxx = 0, u(t, 0) = κ(t), u(t, L) = 0, u x (t, L) = 0, u(0, x) = u 0 (x),

where κ denotes the control input and u 0 is the initial condition. Some full state feedback controls have already been designed in the literature. Let us mention [START_REF] Cerpa | Rapid exponential stabilization for a linear Korteweg-de Vries equation[END_REF] where a Gramianbased method is applied in the case where the control acts on the right endpoint and [START_REF] Tang | Stabilization of Linearized Korteweg-de Vries Systems with Anti-diffusion[END_REF] where the backstepping method is applied with different boundary conditions. With the same approach than the last paper, we have [START_REF] Cerpa | Rapid stabilization for a Korteweg-de Vries equation from the left Dirichlet boundary condition[END_REF] where system (2) is considered.

However, in most cases, we have no access to the full state of the system, and it is more realistic to design an output feedback control, i.e., a feedback law depending only on some partial measurements of the state.

For linear and autonomous finite-dimension systems, a stabilizability and an observability assumptions are sufficient to ensure that the separation principle holds. In other words, if there exists a controller, which asymptotically stabilizes the origin of the system and an observer which converges asymptotically to the state system, the output feedback built from this observer and this state feedback asymptotically stabilizes the origin of the system. In a PDE framework this principle is no longer true and the stability of the closed-loop system is not guaranteed.

The basic question to state the problem is which kind of measurements are we going to consider. The case of a boundary measurement is the most challenging one and the natural choice for the KdV equation ( 2) should be y(t) = u x (t, 0). Unfortunately, the system is not observable with this choice. In fact, if L = 2π, then the stationary solution u(t, x) = 1 -cos(x), satisfies (2) with κ(t) = 0 and in addition u x (t, 2π) = 0 for any time t. The length 2π is not the only one for which observability does not hold (see e.g. [START_REF] Rosier | Exact boundary controllability for the Korteweg-de Vries equation on a bounded domain[END_REF], [START_REF] Cerpa | Control of a Korteweg-de Vries equation: a tutorial[END_REF]).

In this paper we consider the output given by

y(t) = u xx (t, L). (3)
By using this measurement, we build an observer and apply the backstepping method to design an output feedback control which exponentially stabilizes the closed-loop system. The trace u xx (t, L) can be considered itself as a boundary condition or even a control ( [START_REF] Colin | An initial-boundary-value problem for the Korteweg-de Vries equation posed on a finite interval[END_REF], [START_REF] Cerpa | Boundary controllability of the Korteweg-de Vries equation on a bounded domain[END_REF]) and therefore the choice of this output is not an artificial one. This paper is organized as follows. In Section II, we state our main result. Section III is devoted to recall the state feedback control designed in [START_REF] Cerpa | Rapid stabilization for a Korteweg-de Vries equation from the left Dirichlet boundary condition[END_REF]. In IV, we state some regularity results needed to consider the measurement (3) as a continuous function. The observer is built in Section V. In Section VI, the stability of the closed loop controllerobserver system is proven. Finally, Section VII states some conclusions.

II. MAIN RESULT

Based on [START_REF] Krstic | Boundary Control of PDEs: A Course on Backstepping Designs[END_REF] and [START_REF] Smyshlyaev | Backstepping observers for a class of parabolic PDEs[END_REF], we construct the following observer:

     ût + ûx + ûxxx + p 1 (x)[y(t) -ûxx (t, L)] = 0, û(t, 0) = κ(t), û(t, L) = ûx (t, L) = 0, û(0, x) = 0. (4)
Theorem 1. For any λ > 0, there exist a feedback law κ(t) := κ(û(t, x)), a function p 1 = p 1 (x), and a constant C > 0 such that the coupled system (2)-( 3)-( 4) is globally exponentially stable with a decay rate equals to λ, i.e., for any u 0 ∈ H 3 (0, L) we have

u(t, •) H 3 (0,L) + û(t, •) L 2 (0,L) ≤ Ce -λt u 0 H 3 (0,L) (5)

III. FULL STATE FEEDBACK DESIGN

In [START_REF] Cerpa | Rapid stabilization for a Korteweg-de Vries equation from the left Dirichlet boundary condition[END_REF] the following Volterra transformation is introduced

w(x) = Π(u(x)) := u(x) - L x k(x, y)u(y)dy. (6)
The function k is chosen such that the trajectory u = u(t, x), solution of ( 2), with control

κ(t) = L 0 k(0, y)u(t, y)dy, (7)
is mapped into the trajectory w = w(t, x), solution of the linear system

w t + w x + w xxx + λw = 0, w(t, 0) = w(t, L) = w x (t, L) = 0, (8)
which is exponentially stable if λ > 0. In fact, from a Lyapunov approach, we get

d dt L 0 |w(t, x)| 2 dx = -|w x (t, 0)| 2 -2λ L 0 |w(t, x)| 2 dx ≤ -2λ L 0 |w(t, x)| 2 dx, (9)
which gives an exponential decay rate equals to λ for the L 2 -norm of the state w.

The kernel function k = k(x, y) is characterized by:

             k xxx + k yyy + k x + k y = -λk, in T , k(x, L) = 0, in [0, L], k(x, x) = 0, in [0, L], k x (x, x) = λ 3 (L -x), in [0, L], (10)
where

T := {(x, y)/x ∈ [0, L], y ∈ [x, L]}.
The solution of ( 10) exist. This is proved by using the method of successive approximations. Unlikely the case of heat or wave equations, we do not have an explicit solution.

In [START_REF] Cerpa | Rapid stabilization for a Korteweg-de Vries equation from the left Dirichlet boundary condition[END_REF] it is proved that the transformation ( 6) linking ( 2) and ( 8) is invertible, continuous and its inverse is also continuous. Therefore, the exponential decay for w, solution of (8), implies the exponential decay for the solution u controlled by [START_REF] Krstic | Boundary Control of PDEs: A Course on Backstepping Designs[END_REF]. Thus, with this method, the following theorem is proven.

Theorem 2. (State feedback stabilization for KdV ([2])

For any λ > 0, there exist a feedback control law κ = κ(u(t, .)) and C > 0 such that u(t, .) L 2 (0,L) ≤ Ce -λt u 0 L 2 (0,L) [START_REF] Tang | Stabilization of Linearized Korteweg-de Vries Systems with Anti-diffusion[END_REF] for any solution of (2)- [START_REF] Krstic | Boundary Control of PDEs: A Course on Backstepping Designs[END_REF].

IV. REGULARITY RESULT

As we said in the introduction, we consider

y(t) = u xx (t, L) (12)
as a partial measurement of the solution. However, since we have the trace of the second derivative with respect to x of u, we need a regularity stronger than in [START_REF] Cerpa | Rapid stabilization for a Korteweg-de Vries equation from the left Dirichlet boundary condition[END_REF]. Indeed, we ask the output y(t) to be a continuous function. Thus we have the following lemma.

Lemma 1. Let us consider system

     u t + u x + u xxx = 0, u(t, 0) = κ(t), u(t, L) = 0, u x (t, L) = 0, u(0, x) = u 0 (x), (13)
where

u 0 ∈ H 3 (0, L) and κ(t) ∈ H 1 (0, T ). Then u ∈ C([0, T ], H 3 (0, L)) ∩ L 2 (0, T ; H 4 (0, L)) and u xx (•, L) ∈ C([0, T ]).
Proof. This proof is based on [START_REF] Glass | Some exact controllability results for the linear kdv equation and uniform controllability in the zero-dispersion limit[END_REF]. Let us consider the following coordinates transformation

v = u t . ( 14
)
The dynamics of v can be written as follows:

     v t + v xxx + v x = 0, v(t, 0) = κ(t) ∈ L 2 (0, T ), v(t, L) = v x (t, L) = 0, v(0, x) = (-u ′′′ 0 -u ′ 0 ) ∈ L 2 (0, L).
(15) By already known well-posedness results for KdV ( [START_REF] Glass | Some exact controllability results for the linear kdv equation and uniform controllability in the zero-dispersion limit[END_REF]), we get

v ∈ C([0, T ], L 2 (0, L)) ∩ L 2 (0, T ; H 1 (0, L)) ⇒u t ∈ C([0, T ], L 2 (0, L)) ∩ L 2 (0, T ; H 1 (0, L)) ⇒u ∈ C([0, T ], H 3 (0, L)) ∩ L 2 (0, T ; H 4 (0, L)) ⇒u xx ∈ C([0, T ], H 1 (0, L)) ∩ L 2 (0, T ; H 2 (0, L)) ⇒u xx ∈ C([0, T ] × [0, L]) ⇒u xx (•, L) ∈ C([0, T ]) (16)
Thus it concludes the proof of Lemma 1.

V. OBSERVER DESIGN

Based on [START_REF] Krstic | Boundary Control of PDEs: A Course on Backstepping Designs[END_REF], and more precisely on [START_REF] Smyshlyaev | Backstepping observers for a class of parabolic PDEs[END_REF], we can write, for system

     u t + u x + u xxx = 0, u(t, 0) = κ(t), u(t, L) = u x (t, L) = 0, y(t) = u xx (t, L), (17)
the corresponding observer

ût + ûx + ûxxx + p 1 (x)[u xx (t, L) -ûxx (t, L)] = 0, û(t, 0) = κ(t), û(t, L) = ûx (t, L) = 0.
(18) The construction of the observer is based on the finitedimensional design for ẋ = Ax + Bu, y = Cx, which proposes the observer ẋ = Ax + Bu + L(y -C x). If we consider the error e = x -x, then ė = (A -LC)e, and we have to look for a matrix L insuring a good performance. Because of the infinite-dimensional framework we are working in, a matrix L is not enough. Thus we need a function p 1 (x).

In our case, we consider the error ũ := u -û, which satisfies

ũt + ũx + ũxxx -p 1 (x)ũ xx (t, L) = 0, ũ(t, 0) = ũ(t, L) = ũx (t, L) = 0. ( 19
)
Given a positive parameter λ, we look for a transformation Π o defined by

ũ(x) = Π o ( w(x)) = w(x) - L x p(x, y) w(y)dy (20)
such that the trajectory ũ, solution of ( 19) is mapped into the trajectory w = w(t, x), solution of the linear system wt + wx + wxxx + λ w = 0,

w(t, 0) = 0, w(t, L) = 0, wx (t, L) = 0, (21)
which is exponentially stable with a decay rate depending on the value of λ as shown in [START_REF] Rosier | Control and stabilization of the Kortewegde Vries equation: recent progresses[END_REF]. Now, the key step is to find the kernel p = p(x, y) such that ũ(t, x) = Π o ( w(t, x)) satisfies (19). By focusing on (20) and using the Leibniz rules, we get:

• Differentiation along ( 21)

ũt = wt (x) - L x p(x,
-p yy (x, x) -λ) + p(x, L) w(L, t) + (p(x, L) -p 1 (x)) wxx (L, t) -p y (x, L) wx (L, t). ( 26
)
From this equation, we get four conditions:

1. Equation for (x, y) ∈ T :

p yyy (x, y) + p xxx (x, y) + p y (x, y) + p x (x, y) = λp(x, y). ( 27
)
2. First boundary condition on (x, x) for x ∈ [0, L]:

2 d dx p(x, x) + p x (x, x) + p y (x, x) = 0. (28)
3. Second boundary condition on (x, x) for x ∈ [0, L]:

d 2 dx 2 p(x, x) + d dx p x (x, x) + p xx (x, x) -p yy (x, x) -λ = 0. ( 29
)
4. Appropriate choice of p 1 :

p(x, L) = p 1 (x). ( 30
) Recall that T := {(x, y)/x ∈ [0, L], y ∈ [x, L]}.
Moreover, note also that, by setting x = 0 in (20), we get:

p(0, y) = 0, ∀y ∈ [0, L]. (31)
In addition, we have

w(t, L) = ũ(t, L) = wx (t, L) = 0.
Finally, the kernel p satisfies the following PDE:

                 p xxx (x, y) + p yyy (x, y) + p y (x, y) + p x (x, y) = λp(x, y), (x, y) ∈ T , p(x, x) = 0, x ∈ [0, L], p x (x, x) = λ 3 x, x ∈ [0, L], p(0, y) = 0, y ∈ [0, L]. ( 32
)
Let us make the following change of variable:

x = L -y, ȳ = L -x, (33)
and define F (x, ȳ) := p(x, y). Hence:

                 F xxx (x, ȳ) + F ȳ ȳ ȳ (x, ȳ) + F ȳ (x, ȳ) + F x(x, ȳ) = -λF (x, ȳ) (x, ȳ) ∈ T F (x, x) = 0 x ∈ [0, L] F x(x, x) = λ 3 (L -x) x ∈ [0, L] F (x, L) = 0 ȳ ∈ [0, L] (34)
This PDE has already been studied in [START_REF] Cerpa | Rapid stabilization for a Korteweg-de Vries equation from the left Dirichlet boundary condition[END_REF], where no explicit solution has been found, but where the existence of a solution has been proved. Hence, we can conclude that the kernel p := p(x, y) exists. Note that the function Π o defined by ( 20) is linear (by definition) and continuous (because of the existence of p).

VI. STABILITY ANALYSIS OF THE CLOSED LOOP SYSTEM

Instead of dealing directly with the controlled state u and the observer state û, we consider the evolution of the couple (ũ, û) where ũ stands for the error ũ = u -û, as introduced in Section V.

By using the output feedback control

κ(t) = L 0 k(0, y)û(t, y)dy, (35)
the transformation Π defined in ( 6) and its inverse, and the transformation Π o defined in (20) and its inverse, we can see that (ũ, û) are mapped into ( w, ŵ) = (Π -1 o (ũ), Π(û)) solutions of the target system

                     ŵt + ŵx + ŵxxx + λ ŵ = -p 1 (x) - L x k(x, y)p 1 (y)dy wxx (t, L), ŵ(0) = ŵ(L) = ŵx (L) = 0, wt + wx + wxxx + λ w = 0, wx (0) = w(L) = wx (L) = 0. (36)
Note that the parameter λ is the same for the observer and the system itself.

Given that Π, Π o are continuous maps, invertible and their inverse maps are also continuous, the exponential stability of (36) would imply the exponential stability of the closed loop system and therefore the proof of Theorem 1 will be ended.

In order to prove the exponential stability of (36), we use a Lyapunov argument. Let us consider the following function,

V (t) = V 1 (t) + V 2 (t) + V 3 (t), (37)
where

V 1 (t) = A 2 L 0 | ŵ(t, x)| 2 dx (38) V 2 (t) = B 2 L 0 | w(t, x)| 2 dx (39) V 3 (t) = B 2 L 0 | wt (t, x)| 2 dx (40)
with A, B to be chosen later.

Remark 1. We can prove that this Lyapunov function is equivalent to the one obtained by replacing V 3 (t) by

Ṽ3 (t) = B 2 L 0 | wxxx (t, x)| 2 dx (41)
By this, we mean that the exponential decay of one of them implies the exponential decay of the other one. In fact, we can prove that there exist positive constants d 1 , d 2 such that

d 1 (V 2 (t) + Ṽ3 (t)) ≤ V 2 (t) + V 3 (t) ≤ d 2 (V 2 (t) + Ṽ3 (t)).

•

Taking the time derivative of the function V (t), we get after some computations that

V1 (t) =A L 0 ŵt (t, x) ŵ(t, x)dx ≤ (-Aλ + D 2 ) L 0 | ŵ(t, x)| 2 dx + A 2 | wxx (t, L)| 2 = 2 -λ + D 2 A V 1 (t) + A 2 | wxx (t, L)| 2
where D := max x∈[0,L] p 1 (x) -L x k(x, y)p 1 (y)dy . From the same computation as in ( 9), we see that

V2 (t) ≤ -2λV 2 (t).
Moreover, thanks to the regularity H 3 (0, L), the same computation can be applied to z = wt (see the proof of Lemma 1) to obtain V3 (t) ≤ -2λV 3 (t) Thus, we get:

V (t) ≤2 -λ + D 2 A V 1 (t) + A 2 | wxx (t, L)| 2 -2λV 2 (t) -2λV 3 (t).
We need to find an upper bound for | wxx (t, L)| 2 . We multiply wt + wx + wxxx + λ w = 0,

w(0) = w(L) = wx (L) = 0, (42)
by x wxx and after some computations we get

| wxx (t, L)| 2 ≤ 1 L + L wxx 2 L 2 (0,L) + 2λ + 1 L wx 2 L 2 (0,L) + 1 L wt 2 L 2 (0,L) (43)
and finally the existence of a, b > 0 such that

| wxx (t, L)| 2 ≤ a w 2 L 2 (0,L) + b wt 2 L 2 (0,L) (44)
Remark 2. Here, we have used that the norm f H 3 (0,L) and the norm f L 2 (0,L) + f xxx L 2 (0,L) are equivalent. See also Remark 1.

•

We use the latter inequality to write:

V (t) ≤2 -λ + D 2 A V 1 (t) + 2a A 2 B V 2 (t) + 2b A 2 B V 3 (t) -2λV 2 (t) -2λV 3 (t).
Therefore,

V (t) ≤2 -λ + D 2 A V 1 (t) + 2 -λ + aA 2 B V 2 (t) 2 -λ + bA 2 B V 3 (t).
In this way, by tuning A, B large enough, we get for any ε > 0 that

V (t) ≤ 2 -λ + ε V (t),
which gives an exponential stability with decay rate as close to λ as we want. The rapid stabilization is achieved because the parameter λ can be chosen as large as desired.

It concludes the proof of the stability of the closed loop system with the output feedback control law depending on a boundary measurement of the state.

VII. CONCLUSIONS

In this paper we have designed an output feedback controller for the linear Korteweg-de Vries equation posed on a bounded interval. This boundary feedback control acts on the Dirichlet boundary condition on the left endpoint and exponentially stabilizes the system by using the boundary measurement y(t) = u xx (t, L). Because of this choice of y(t), we have to work on a more regular framework, given by H 3 (0, L) as the state space. The backstepping method is applied together a classical observer in order to build the control.

This work is the first one addressing the output feedback problem for the Korteweg-de Vries equation and opens several possible extensions. Different location for the control or different boundary conditions can be considered as in [START_REF] Cerpa | Rapid stabilization for a Korteweg-de Vries equation from the left Dirichlet boundary condition[END_REF], [START_REF] Tang | Stabilization of Linearized Korteweg-de Vries Systems with Anti-diffusion[END_REF]. The choice of the measurement is particularly interesting. For instance, it would be nice to deal with the collocated case, which could be harder than the noncollocated problem considered in this paper. Of course, the nonlinear case also appears as a natural next step in the study of this control system. In that context, even performing simulations of the closed loop system may be an interesting challenge.

*This work has been partially supported by FONDECYT grant 1140741, MathAmSud COSIP, and CONICYT grant ACT-1106


Zhibin Zhou
email: zhibin.zhou@ecole-navale.fr

Franck Scuiller

Jean Frédéric Charpentier

Mohamed Benbouzid
email: mohamed.benbouzid@univ-brest.fr

Tianhao Tang
email: thtang@shmtu.edu.cn

An Up-to-Date Review of Large Marine Tidal Current Turbine Technologies

Keywords: Tidal current turbine, large machine, review

Owning to the predictability of tidal current resources, marine tidal current energy is considered to be a reliable and promising renewable power source for coastal areas or some remote islands. During the last 10 years, various original horizontal axis and vertical axis marine current turbines (MCT) have been developed around the world. Although various projects have been reported in the state-ofthe-art research papers in recent years, many of these projects were only at the design stage when the papers were published. In fact, some projects do not have any further developments during the several years after the first reporting. In this paper, up-todate information about large tidal turbine projects over 500 kW is focused. The newest achievements of these large tidal current turbine technologies are presented. These technologies represent the industrial solutions for several pre-commercial MCT farm projects in the coming years. This paper provides a useful background for researchers in the marine turbine energy domain.

I. INTRODUCTION

One of the main advantages of marine current energy is related to the predictability of the resource. The astronomic nature of tides is driven by the gravitational interaction of the Earth-Moon-Sun system and makes marine tidal currents highly predictable with 98% accuracy for decades [START_REF] Benelghali | A simulation model for the evaluation of the electrical power potential harnessed by a marine current turbine[END_REF]. There are basically two ways of generating electricity from marine tidal energies: either by building a tidal barrage across an estuary or a bay, or by extracting energy from free flowing tidal currents. The main drawback of the solution is that large barrage system would change the hydrology and may have negative impacts on the local ecosystem [START_REF] Pelc | Renewable energy from the ocean[END_REF]. Therefore during the last few decades, developers have shifted towards technologies that capture the kinetic energy from tidal-driven marine currents. The exploitable marine current energy with present technologies is estimated about 75 GW in the world and 11 GW in Europe [START_REF] Boye | Rapport de la mission d'étude sur les énergies marines renouvelables[END_REF].

In fact, various original horizontal axis and vertical axis marine current turbines (MCT) have been developed around the world in recent years [START_REF] Rourke | Marine current energy devices: Current status and possible future applications in Ireland[END_REF][START_REF] Bahaj | Generating electricity from the oceans[END_REF][START_REF] Zhang | State of the art in tidal current energy extracting[END_REF]. The majority of MCT devices are horizontal axis turbines with rotation axis parallel to the This work is supported by Brest Métropole Océane (BMO) and the Shanghai Maritime University. current flow direction. The main disadvantages associated with vertical axis turbines are relative low self-starting capability, high torque fluctuations and generally lower efficiency than horizontal axis turbine design. Currently, only horizontal axis MCTs appear to be the most technologically and economically solution for large-scale marine current turbines with power capacity over 500 kW.

Although various turbine projects have been reported in some state-of-the-art research papers, many of these projects were only at the design stage when the papers were published [START_REF] Rourke | Marine current energy devices: Current status and possible future applications in Ireland[END_REF], [START_REF] Benelghali | Marine tidal current electric power generation technology: State of the art and current status[END_REF]. However, some projects are abandoned or never have been built, for instance, the Lunar Energy Tidal turbine. And some projects do not have further developments during the several years after the first announcement, such as Atlantis Resource Nereus and Solon turbines. Therefore, an up-to-date review of industrialized turbine technologies is necessary.

In this paper, up-to-date information and the newest achievements of industrialized large tidal current turbine technologies (over 500 kW) will be focused. In Section II, the pilot pre-commercial MCT farms information and the industrialized large turbine technologies are presented. The conclusion and perspective are then given In Section III.

II. LARGE MARINE CURRENT TURBINE TECHNOLOGIES

Several horizontal axis turbine technologies are developed more than one or two generations and have been chosen by the industrial communities to realize pilot demonstrative MCT farms before the final commercial stage. These projects illustrate the up-to-date developments of large MCT technologies which will provide electricity to coastal or island areas in the coming years. Table I summarizes the main information about some of these pilot MCT farm projects and their planned/estimated operational dates. From this table, it can be seen that most of these turbine technologies have attended megawatt-level power capacity.

A. OpenHydro Turbine Technology

OpenHydro is an open-center turbine technology; a 250 kW prototype was installed and tested at European Marine Energy Center (EMEC) off Orkney islands in Scotland and was connected to the UK national grid in 2008. This turbine

B. HS1000 and AR1000 Turbine Technologies

The 1MW pre-commercial turbine HS1000 (Fig. 2) was tested by Andritz Hydro Hammerfest (original Hammerfest Strøm) at EMEC tidal test site at the end of 2011. It started delivering energy to the grid in 2012. The HS1000 turbine is based on the technology of a smaller prototype HS300 (300 kW, reported as E-Tide turbine project in [START_REF] Zhang | State of the art in tidal current energy extracting[END_REF][START_REF] Benelghali | Marine tidal current electric power generation technology: State of the art and current status[END_REF]) which was installed in Norway and connected to the public grid in 2004. This megawatt-level HS1000 turbine technology is planed to be used in a 10 MW commercial array in the Sound of Islay on the west coast of Scotland; the tidal resource and seabed surveys are completed [13]. This technology is also reported to be chosen in the first phase of the MeyGen tidal current project in Inner Sound of the Pentland Firth [9].

Another candidate for the MeyGen project is the AR1000 turbine (Fig. 3) technology developed by Atlantis Resources Corporation [14]. The Atlantis AR1000 turbine features fixed pitch configuration and is rated at 1 MW at 2.65 m/s current velocity. The AR1000 can be rotated in the slack period between tides using a yaw drive for optimally facing the tides. The first AR1000 was successfully deployed and tested at the EMEC facility during the summer of 2011. The AR1000 system is also scheduled to be installed on Daishan demonstration site in China during 2014. A larger turbine called AR1500 (1.5 MW at 3.0 m/s) is under development for future installation in the Pentland Firth in Scotland and the Bay of Fundy in Canada.

C. SeaGen S Turbine Technology

The SeaGen S turbine (showed in Fig. 4) developed by Marine Current Turbine Ltd. (owned by Siemens since 2012) is a well-known system that can generate electricity from marine current energy. This technology has a twin axial-flow turbine supported on a structure with the ability to raise the moving components out of the water for maintenance [16]. It is the world first grid-connected megawatt-level MCT system. The 1.2 MW SeaGen S system (2×600 kW) was installed in Strangford Lough in Northern Ireland in 2008 and has generated 8 GWh electricity since the installation. During a strong spring tide, the SeaGen S system can deliver more than 20 MWh in a single day. It is reported that up-scaled SeaGen S systems with 20 meter rotor and 2 MW power rating will be installed in two commercial arrays in UK waters (as listed in Table 1) from 2015. The Voith Hydro turbine adopts robust design: it uses fixed-pitch blades and permanent magnet synchronous generator to achieve compact structure. Specially designed symmetric blade profile can be operated for bidirectional tidal currents avoiding pitch and yaw requirements. The nacelle system uses seawater lubrication technology to avoid grease and seals for the bearings. The electrical machine is protected from corrosion and marine fouling by proven coatings and sacrificial anodes. The low-maintenance characteristics makes the estimated system lifetime up to 20 years [18]. The 1 MW turbine has a rotor diameter of 16 m and reaches the rated power at a current velocity of 2.9 m/s.

D. Voith Hydro Turbine Technology

E. Sabella Turbine Technology

One of the potential MCT applications is to provide electricity for remote islands not connected to the continental electric grid. As example, in the near future, several MCTs are planed to be installed in the Fromveur passage (near Brest) in France to supply part of the load demand of Ouessant Island. For the Fromveur project, the Sabella tidal turbine technology will be used [19]. The Sabella D10 turbines have a rotor diameter of 10 m and a power capacity of 0.5~1.1 MW for 3.0~4.0 m/s current velocities. Figure 6

F. Alstom Tidal Turbine Technologies

In 2010, Alstom announced a 1 MW turbine project named Beluga 9 for deployment in Canada's Bay of Fundy in 2012 [21]. Beluga 9 was a joint project with Canadian company Clean Current. It was intended for very powerful currents (up to 4.5 m/s) and for sites with depths of 30 meters or more. However, with the termination of the licensing agreements with Clean Current at the end of 2012, there is no further news about this megawatt turbine project. Clean Current has decided to focus on river turbines and shallower water tidal turbines (< 20 meters). While Alstom's intention is to pursue the larger utility scale opportunities [22].

In 2013, Alstom successfully installed a 1 MW tidal current turbine at EMEC tidal test site after the acquisition of Tidal Generation Limited (TGL) [23][24]. Figure 7 shows the Alstom's 1 MW tidal turbine system. This tidal turbine weighs 150 tonnes and has three pitchable blades, a rotor diameter of 18 m and a nacelle length of 22 m. A standard drive train and power electronics device is inside the nacelle. This turbine system can be installed in a water depth between 35 and 80 meters. It can generate electricity with a tidal current between 1 m/s and 3.4 m/s, and reaches the rated power for a tidal current speed of 2.7 m/s. This 1 MW Alstom tidal turbine is under extensive testing and analysis in different operational conditions off Orkney islands throughout 2013 over an 18 month period. Although the Alstom's 1 MW tidal current turbine seems similar to a standard wind turbine system (with three pitchable blades and gearbox-integrated drive train), this tidal turbine have some special characteristics. Its buoyancy enables the turbine to be easily installed and retrieved by using small vessels, thus reducing the installation and maintenance costs. It has an intelligent nacelle which can be rotated by the water thrusts to orientate the turbine system to face the tide direction; that enables the system to harness energy efficiently from both ebb and flood tides

III. CONCLUSIONS AND PERSPECTIVES

The turbines presented above represent the newest achievements in the megawatt-level marine current turbine technologies. The comment point is that they are all horizontal axis turbine. And it should be noted that more than half of these turbine technologies adopt the fixed pitch solution for the blades.

In addition to these industrialized large turbine technologies, there are other megawatt turbine projects under development. For example, the TidalStream Triton technology mounts several turbines on a semi-submersible buoying platform to achieve high power harness capacity. It should be noted that in FORCE, there are also several megawatt-level current turbine plans which are strongly related to the turbine technologies presented in this paper. High tidal current speed (up to 5.5m/s) and Nova Scotia's feed-in-tariff make FORCE one of the most attractive and economic tidal sites in the world [26][27]. Siemens and Bluewater are jointly developing a 2 MW floating tidal current turbine, called SeaGen F, to be installed in Bay of Fundy and it will supply electricity for 1,800 Nova Scotia households. OpenHydro, the DCNS tidal subsidiary, will proceed with plans for a grid-connected 4 MW tidal array to be installed at FORCE in later 2015. This array will consist of two 2 MW turbines with 16 m rotor diameters.

  technology has been chosen by the French companies EDF and DCNS to build a demonstrative MCT farm off the coast of Paimpol-Bréhat in Brittany, France. The first 500 kW OpenHydro turbine (Fig.1) was tested in September 2011 near Brest. This 850 tonnes turbine has a diameter of 16 m and is supposed to be installed at a depth of 35 meters. This technology uses multi-blades fixed between the open-center rim and the outside shell in a rim-driven configuration. The permanent synchronous generator is integrated into the outside rim shell. The planed MCT farm with 4 turbines is reported to be in operation in 2014[10]. However, some delays in the final farm operation could still be envisaged.

Fig. 1 .

 1 Fig. 1. DCNS-OpenHydro turbine © [11].

Fig. 2 .

 2 Fig. 2. Andritz Hydro Hammerfest HS1000 turbine © [12].

French

  energy company GDF Suez has plans to install pilot tidal energy farms at Raz Blanchard off the coast of Lower Normandy and in the Fromveur passage off the coast of Finistère in Brittany. These two sites represent 80% of the marine current energy potential in France [17]. For the Raz Blanchard project, GDF Suez has recently confirmed to use Voith Hydro HyTide turbine and Alstom tidal turbine technology. The Voith Hydro turbine system (shown in Fig. 5) is developed by German hydropower equipment maker Voith Company. The first test turbine of 110 kW has been in operation near the South Korean island of Jindo since 2011. The up-scaled version of 1 MW turbine is now installed and tested at EMEC tidal test site.

  illustrates the perspective Sabella turbine farm. The Sabella D10 is based on the first French marine current turbine Sabella D03 (3 m rotor diameter) which was tested at Bénodet estuary near Brest in 2008. The prototype D10 turbine is now completing the construction and is scheduled to be installed in the Fromveur passage at the end of 2014. Larger turbines D12 and D15 with power capacities of 1~2 MW are under design for future turbine farm applications [19-20].

Fig. 3 .

 3 Fig. 3. Atlantis AR1000 turbine © [14].

Fig. 4 .

 4 Fig. 4. SeaGen S turbine © [15].

Fig. 5 .

 5 Fig. 5. Voith Hydro turbine © [18].

Fig. 6 .

 6 Fig. 6. Sabella turbine farm illustration © [20].

Fig. 7 .

 7 Fig. 7. Alstom tidal turbine © [23].

  The 1/23rd scale and 1/10th scale turbines (Triton T6 and Triton T3) were successfully tested in 2009 and 2011 respectively [25]. The first full-scale version Triton T36 system (with 2.5 MW power rating) is planned to be built at the Fundy Ocean Research Center for Energy (FORCE) tidal test site in Nova Scotia, Canada.

TABLE I PILOT

 I MCT FARMS IN THE COMING YEARS [3], [8-10].

	Companies	Location	Turbine Name	Pitchable Blades	Turbine Number	Total Capacity (MW)	Operational Year
	DCNS, EDF	Paimpol-Bréhat	OpenHydro	No	4	2	2014/2015
	MeyGen	Pentland Firth (Scotland)	HS1000 or AR1000	Yes No	6	6	2015/2016
	MCT, Siemens	Kyle Rhea (Scotland)	SeaGen S	Yes	4	8	2015
		Anglesey (Wales)	SeaGen S	Yes	5	10	2015
	Andritz Hydro Hammerfest	Sound of Islay (Scotland)	HS1000	Yes	10	10	N/A
	GDF Suez, Eole	Raz Blanchard	Voith Hytide	No	3~6	3~12	2016
	Generation	Fromveur	Sabella	No	N/A	N/A	2016


Sofiane Djebarri
email: sofiane.djebarri@ecole-navale.fr

Jean Frederic Charpentier

Franck Scuiller

Mohamed Benbouzid
email: mohamed.benbouzid@univ-brest.fr

Comparison of Direct-Drive PM Generators for Tidal Turbines

Keywords: Marine Current turbine, POD, Rim-Driven, PM Generators

The aim of the paper is to compare the design of direct-drive permanent-magnet (PM) generators associated with horizontal axis tidal turbines. The turbine/generator couplings are examined. These turbine/generator couplings consist of a POD and Rim-driven assembly. In a Rim-Driven association the electrical generator active parts are inserted in a duct surrounding the blades. For POD generator insertion, the electrical machine is placed in a nacelle located on the turbine axis. To achieve the generators sizing, a design rated point related to an industrial MCT is defined. The used design models include an electromagnetic model which is linked to a thermal model in an optimization procedure that goals to minimize the active parts cost. Firstly, a single rotor/single stator PM axial flux generator and a radial flux PM generator are designed for a rimdriven MCT specification. For these generators sizing, a comparison of the machines active parts and the machines geometrical dimensions are carried out. Secondly, radial flux PM generators are designed for Rim-Driven and POD assembly and a comparison is performed for this study case. Finally, the influence of the POD diameter on the generator electromagnetic design is studied. It shows that the active parts costs are minimized, when the generator diameter is around 1/3 of the turbine diameter for the considered specifications. These performed comparisons between the considered study cases aims to help designers in their technologies choices.

I. INTRODUCTION

Direct driven PM generators are relevant solution for marine and offshore energy applications because they allow minimizing the maintenance operations on the systems. In direct driven systems, generator is operated at low speed (typically lower than 50rpm). In one hand, the elimination of the gearbox leads to an improvement of efficiency and a minimization of maintenance [START_REF] Semken | Direct-drive permanent magnet generators for high-power wind turbines: Benefits and limiting factors[END_REF]. In the other hand the low speed generator is characterized by a high torque. This high torque level is related to high active parts volume (iron, magnets, copper).

Several hydrodynamic principles are used in the major current projects [START_REF] Block | Tidal power: an update[END_REF]: horizontal [START_REF] Drouen | Investigations on the performances of the electrical generator of a rim-driven marine current turbine[END_REF][4][START_REF] Gilson | Alstom's tidal in stream energy converter demonstration project[END_REF] or vertical axis turbines or with oscillating hydrofoils [START_REF] Benelghali | Marine tidal current electric power generation technology: State of the art and current status[END_REF]. These solutions can be associated with geared systems or direct driven systems [START_REF] Benelghali | Marine tidal current electric power generation technology: State of the art and current status[END_REF]. The developed solutions are all at a TRL level between 5 and 7 and reference solutions have not already been determined at the industrial level. This is why it is necessary to investigate the pro and the cons of the possible technological solutions.

The presented work presents elements of comparison for several direct driven generators for horizontal axis turbine. The possibility to locate the generator in a duct which surrounds the turbine (rim driven system) like in the Open Hydro project or in [START_REF] Drouen | Investigations on the performances of the electrical generator of a rim-driven marine current turbine[END_REF] (Fig. 1) is compared to the more classical solution where generator is located in a pod located in the turbine axis (fig. 2) like in the Atlantis Project or the abandoned ALSTOM Beluga 9 project [START_REF] Gilson | Alstom's tidal in stream energy converter demonstration project[END_REF], [7].

In the first case, the internal radius of the generator is equal to the external radius of the turbine blades. In the second case, the external radius of the generator is limited to avoid too strong disturbance in the flow in the turbine disk.

II. DESIGN METHODOLOGY

A. Used Methods and Specifications

The aim of the study is to compare and evaluate several topologies of tidal current generator in the EM structure point of view (axial and radial flux ) and in terms of assembly point of view (rim-driven or pod assembly). The design methodology couples an EM model and a thermal model. These two models are inserted in an optimization process which aims to minimize the active parts cost. Constraints are considered in the optimization process. They are related to maximal temperature, minimal efficiency and geometry (internal and external radius and slot shape). The EM model is derived from a 2D model proposed in [START_REF] Zhu | Instantaneous magnetic field distribution in brushless permanent magnet dc motors, Parts I to IV[END_REF] and solves the governing field equations by separating the polar variables. This EM model allows the calculation of the dimensions and performance of the generator from EM material characteristics (saturation level of soft materials, magnet magnetization) and operating point specifications. The thermal model, is mainly similar to the model presented in [START_REF] Djebarri | Rough design of a double-stator axial flux permanent magnet generator for a rim-driven marine current turbine[END_REF], and allows checking that the temperature does not exceed a maximal limit. Common data are used. They relate to slot filling factor, to winding topology (diametral windings), to ambient temperatures and to active material characteristics. Common site and turbine specifications are also used. A specific tidal speed is determined from tidal current speed statistics on a typical tidal spot. A typical optimal turbine power coefficient corresponding to an optimal Tip Speed Ratio (TSR) is fixed (these data are related to turbine hydrodynamic performance). With the knowledge of the turbine diameter, these data allows defining a characteristic rotating speed and corresponding power which corresponds to the common design operating point which is used in all the studied cases for the generator designs. This common data and specifications allow making relevant comparison between the studied design options. These data are summarized in table I and table II. The comparison results are related to the global active part cost and mass. A scheme of the association of the models and optimization method is given at Fig. 3.

B. Studied Magnetic Structures

The axial flux machine under consideration is a mono stator mono rotor structure similar to the one depicted in fig. 4 (from [START_REF] Aydin | Axial flux permanent magnet disc machines: A review[END_REF]). This structure comprises a stator disk which supports windings and a rotor disk supporting magnets. Of course, the active conductors are oriented along radius direction and the magnet flux is oriented in the axial direction. The second structure is a more classical radial flux structure like the one presented in fig. 5. In this structure the active conductors are oriented along the axial direction and the magnet flux along the radial direction. For high diameter generators a hollow center can be used. These two kinds of structure are associated with the studied turbine in pod type assembly or in a rim driven assembly. For instance, Fig. 6 shows a scheme of a radial flux machine in a rim driven turbine.

C. Site and Turbine Characteristics

The used design operating point relates to a 1MW turbine located in a high tidal current spot. The operating point corresponds to a tidal current speed of 3.5m/s. The studied turbine is characterized by an optimal power coefficient of 0.41 corresponding to a TSR (Tip Speed Ratio) of 6. This value corresponds to a 34rpm rotating speed and a 1MW power for a 12m turbine diameter. These data are summarized in table III. In this section, a radial flux PM (RFM) machine and an axial flux PM (AFM) machine are designed for a Rim-Driven assembly considering the methodology and common specifications described in the previous sections.

Table IV shows the obtained results regarding the geometrical dimensions and the performances. 7b show a sketch of a few poles of the two structures (axial view for the RFM and radial view for the AFPM). Figure 7a shows the overall dimensions of the RFM. These 3 figures show the main dimensions of the two structures.

Figures 8a and8b give the calculated cost and mass for the active parts of the 2 Rim-Driven structures.

As shown in Fig. 8, the cost and mass of the active parts of the two structures are roughly of the same order. The main dimensions of the two structures are similar (Table IV). In the two cases, the structures are characterized by a very small axial length which leads to a very high volume of end windings. The end windings copper represents more than 60% of the total copper volume. The constraint on efficiency is saturated and the machine remains relatively cold (the maximal temperature is not reached). This thermal behavior is related to the machine immersion which contribute to a very good cooling. This first part of the study shows that AFM and RFM designs are roughly similar. This is why only Radial Flux Machines are considered in the last part of the paper. On the first hand, results show that the cost of the POD generator is 20% lower than the one of the rim driven structure. This point stigmatizes a better use of magnet and copper relative to a higher value of active length. On the other hand, the mass of the pod generator is higher than the one of the rim driven structure. This higher mass is mainly related to iron parts.

In a second step, optimal designs are calculated with different fixed values of maximal radius of the generator to evaluate the influence of this radius on the design. This value is varying between 2m and 12m. It can be noted that the 12m value corresponds roughly to a rim-driven configuration and that the values between 7 and 11m seem not very realistic in terms of hydrodynamic behavior. The results in terms of cost, mass, active length and pole numbers are presented in Figs. 10a to 10d. Figure 10a shows that an optimal value of the external radius which is around 1/3 of the studied turbine diameter leads to a minimal cost of active parts.

Figure 10b shows that upper a given value around 5m the mass remains roughly constant. This fact can be explained with the reduction of axial length which leads to higher end windings volume and higher flux leakage phenomenon. A too strong external radius constraint (external diameter lower than 4m) leads to a strong rise of the cost and mass Indeed if we consider a given force density  t (N/m 2 ) in the gap, the torque is directly proportional to the air gap surface and the diameter (

1 2 t T SD 
).That means that if end effects and flux leakage effects are neglected the size of the machine is in first order inversely proportional to the diameter as underlined in [START_REF] Drouen | Investigations on the performances of the electrical generator of a rim-driven marine current turbine[END_REF]. This is why a too strong constraint on internal diameter leads to an oversizing of the generator dimensions. In figure 10-d it can be noted that the number of pole pairs increase with the diameter. This phenomenon is related to the fulfillment of a geometrical constraint on the shape of the slots (the ratio between slot height and slot width is limited to a maximal value (Table I).

V. CONCLUSION

The presented study focuses on the design of direct driven permanent magnet generators for tidal current turbine. The aim of this study is to bring out comparison elements between several electromagnetic structures and generator/turbine association. Two electromagnetic topologies are considered: a single rotor/single stator axial flux PM generator and a more classical radial flux generator. Two kind of turbine/generator association are studied: a rim driven system where the generator is included in a duct surrounding the turbine blades and a more classical solution where the generator is located in a pod connected to the turbine shaft. Generator dimensions and performances are calculated for these solutions. These designs correspond to a common rated operating point related to a 12m diameter turbine located in a tidal spot with high current values. The design method associates an EM model and a thermal model in an optimization process which aims to minimize the cost of active parts with respect to geometrical, thermal and efficiency constraints. The study shows that the single rotor single stator axial flux machine cost and mass are of the same order than those obtained for a radial flux machine. Considering the turbine/machine association a POD structure where the radius of the generator is limited to 1/3 of the structure is heavier than the Rim Driven solution. However this POD solution will lead to a smaller global cost of the active parts. The study shows that an optimal value of the external radius of the generator can be found which corresponds to mass or cost minimization for each given specification set (rated operating point and given turbine).

ForFig. 1 .Fig. 2 .

 12 Fig. 1. Rim-Driven experimental turbine IRENav (a), Open Hydro RimDriven (b) [4].

Fig. 3 .

 3 Fig. 3. Models and methods association.

Fig. 4 .

 4 Fig. 4. Basic sketch of a mono stator-mono rotor axial flux machine [13].

Fig. 5 .

 5 Fig. 5. Radial flux machine sketch.

Fig. 6 .

 6 Fig. 6. Basic sketch of a radial flux rim driven generator.

Fig. 8 .

 8 Fig. 7. View of the designed structures

Fig. 9 .

 9 Fig. 9. Comparison of RFM POD and RD structures: (a) active part cost, (b) active part mass.

Fig. 10 .

 10 Fig. 10. POD external diameter influence, (a) active part cost, (b) active part mass, (c) active lenght, (d) pole numbers.

TABLE I .

 I COMMON DESIGN DATA FOR THE EM STRUCTURES.

	Iron laminations FeSi (3% to 3.5% of Si) [12]
	Operating electrical frequency felec	50 to 400	Hz
	Saturation level	Bsat	1.48	T
	Maximal field density	-	2.2	T
	Lamination thickness	-	0.1	mm
	Specific losses (50Hz)	PFe0	5	W/kg
	Relative permeability	-	5000	-
	Density	ρFeSi	7700	kg/m 3
	cost [12]	Clam	1	$/kg
		Copper (20°C)
	Electric resistivity	ρCu	1.6779  10 -8 Ω.m
	density		8960	Kg/m 3
	cost [11]	Ccopper	7.8	$/kg
				Maximum tooth height to tooth width ratio	max 10	-
				Minimum tooth height to tooth width ratio	min	4	-
				Relative air gap thickness related to air gap diameter (D)	k	2‰	-
				Air gap	hg	k.D	mm
				Angular shift between EMF and currents (MTPA sinusoidal control of stator currents)	ψ	0	rad
				Phase number	m	3	-
				Number of slot per pole and per phase	Spp	1	-
				End winding length [10]	Lew	πD/2p m
				Slot filling factor	kf	0.65	-
				Relative magnet width (relative to the pole width) βm	0.66	-
				Maximal copper temperature	Tmax	100	°C
				Minimal electrical efficiency	ηélec	0.90	-
				Sea water temperature	Twater 30	°C
				TABLE II.	ACTIVE PARTS CHARACTERISTICS AND COSTS.
				Permanent magnets (NdFeB) [11]
				Magnetization	Br	1.22	T
				Coercitive field	Hcj	1208	kA/m
				Maximal temperature	-	100	°C
				density	7400	kg/m 3
				cost [11]	Cmagnets 115	$/kg

TABLE III .

 III OPERATING POINT DATA.

	III. RADIAL AND AXIAL FLUX STRUCTURE COMPARISON FOR A
	RIM DRIVEN TURBINE
	Rated value of the water current Vr	3,5	m/s
	Turbine rated power	Pr	1000	kW
	Turbine rated rotating speed	Ωr 34	tr/min
	Rated torque	Tr	280860 Nm
	Turbine diameter	D0 12	m

TABLE IV

 IV

	.	MAIN DIMENSIONS AND EM PERFORMANCE OF THE TWO
			STRUCTURES.
			AFM	RFM
	External radius of the generator	6.0739	6.094	m
	Average radius of the active parts (for AFM only)	6.0369	-	m
	internal radius of the stator (for RFM only)	-	6.0412	m
	Internal radius=turbine external radius generator	6	6	m
	Active lenght		7.39	7.75	cm
	Number of pole pairs	232	238	-
	Electric efficiency (ηelec)	0.96	0.96	-
	Iron losses percentage (Pfe)	22.17	22.87	%
	Joule losses percentage (PJ)	77.82	77.12	%
	Conductor temperature	maximal	50	53.7	°C
	Ratio of the end windings
	copper volume on the total	63.45	61.88	%
	copper volume.
	Figures 7c and


Laurent Schwaederlé

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

Cathy Rond

Alix Gicquel

Molecular dynamics calculations of CH 3 sticking coefficient onto diamond surfaces

Keywords: molecular dynamics, deposition, sticking, diamond, methyl 1

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

In power electronics, the semiconductor switching device is a key component because it is a part of all electrical power conversion and control systems [START_REF] Zhe Chen | A review of the state of the art of power electronics for wind turbines[END_REF]. Whereas silicon-based devices still dominate the power electronics applications, Si technology shows difficulty in meeting the demand for some high-power utility applications as a result of limitations in its intrinsic material properties [START_REF] Achard | Thick boron doped diamond single crystals for high power electronics[END_REF]. These limitations as well as the need of weight and volume reduction of power electronics systems require new materials. The wide band gap materials represent the best candidates since they possess better properties than silicon such as breakdown voltage, lower conduction resistance, ... Among them, diamond offers the best combination of properties.

Today, microwave (MW) plasma assisted chemical vapor deposition (MW-PACVD) reactors are the most efficient reactors to produce thick (few 100's µm) high quality / high purity crystals [START_REF] Achard | The control of growth parameters in the synthesis of high-quality single crystalline diamond by CVD[END_REF][START_REF] Achard | High quality MPACVD diamond single crystal growth: high microwave power density regime[END_REF]. However, high power / high pressure conditions are necessary to dissociate efficiently the H 2 -(1-10%) CH 4 gas mixture. Because of the high gas temperatures reached in the plasma bulk (2800 -3500 K) [START_REF] Hassouni | Modelling of diamond deposition microwave cavity generated plasmas[END_REF], the efficient generation of reactive species, in particular atomic hydrogen, a key specie for growing diamond under these conditions, relies on the thermal dissociation of molecular hydrogen rather than on its electron impact dissociation.

Despite several works over the years, the identity of the precursor species for diamond growth is still under discussion. On the one hand, based on either theoretical approach [START_REF] Larsson | Nearest-neighbor influence on hydrocarbon adsorption on diamond (111) studied by ab initio calculations[END_REF][START_REF] Frenklach | Growth mechanism of vapordeposited diamond[END_REF][START_REF] Chang | Minimumenergy paths for elementary reactions in low-pressure diamond-film formation[END_REF][START_REF] Peploski | Molecular dynamics studies of elementary surface reactions of acetylene and ethynyl radical in low-pressure diamond-film formation[END_REF][START_REF] Rau | Modeling of diamond growth from a microwave plasma: C2h as growth species[END_REF] or measurements of the gas com-position near the substrate [START_REF] Dek | Ethynyl (c2h): A major player in the chemical vapor deposition of diamond[END_REF][START_REF] Cappelli | An investigation of diamond film deposition in a premixed oxyacetylene flame[END_REF], it has been suggested that either C 2 H 2 or C 2 H could play a significant role on diamond growth mechanisms. Some molecular dynamics simulations have been carried out in order to reproduce the surface processes for these two species [START_REF] Neyts | Molecular dynamics simulations for the growth of diamondlike carbon films from low kinetic energy species[END_REF][START_REF] Trskelin | Molecular dynamics simulations of c2, C2H, C2H2, C2H3, C2H4, C2H5, and C2H6 bombardment of diamond (1 1 1) surfaces[END_REF]. However, in microwave plasma reactors operating in a moderate pressure range (10 -300 mbar) and at input microwave powers below 5 kW, methyl radical has been recognized as the main carbon species responsible for diamond growth [START_REF] Lombardi | Study of an h 2/CH4 moderate pressure microwave plasma used for diamond deposition: modelling and IR tuneable diode laser diagnostic[END_REF][START_REF] Hassouni | Modelling of diamond deposition microwave cavity generated plasmas[END_REF]. On the other hand, Harris and Goodwin have proposed a growth mechanism on the (100) diamond surface, which is the preferred face to growth high quality single crystals [START_REF] Goodwin | Scaling laws for diamond chemicalvapor deposition. i. diamond surface chemistry[END_REF][START_REF] Goodwin | Scaling laws for diamond chemicalvapor deposition. ii. atomic hydrogen transport[END_REF]. The growth model assumes two main species: the atomic hydrogen H and the methyl radical CH 3 that both play key roles. In the model detailed below, the star identifies unsaturated radicals or molecules and Cd means a carbon atom that will be incorporated in the diamond lattice:

             C d -H + H → C * d + H 2 (a) C * d + H → C d -H (b)                              C * d + CH 3 → C d -CH 3 (c) C d -CH 3 + H → C d -CH * 2 + H 2 (d) C d -CH * 2 + H → C d -C d -H + H 2 (e)
The H atoms 1) stabilize the surface by saturating the dangling bonds, 2) create growth sites, 3) etch the graphitic phase and 4) dehydrogenate the adsorbed radical, while CH 3 is the species responsible for diamond film growth.

Due to the highly dissociated H 2 -CH 4 plasma, the diamond surface dangling bonds are all saturated by hydrogen atoms giving rise to a surface with sites presenting only sp 3 for C-atom configuration. Since the sticking of CH 3 on the surface can only take place in a free radical site (dangling bond), the first step must be the creation of radical site that occurs by hydrogen abstraction when hydrogen atoms impact the surface forming H 2 molecules that leave the surface (Eley-Rideal mechanism) (a). The most probable fate of this newly created radical site is to be quickly hydrogenated again due to the high atomic hydrogen density and mobility (b). The number of unsaturated dangling bonds on the diamond surface, at any one time, is a function of the equilibrium between these two processes (H abstraction and H addition reactions), which are both sensitive functions of the substrate and gas temperatures. Although the probability is much lower than that of H-atoms, CH 3 can react with the surface dangling bonds and stick to them (c). Subsequent dehydrogenation reactions (d) -(e) lead to incorporation of C atoms into the diamond lattice with a sp 3 configuration.

The fraction of dangling bonds (non-hydrogenated surface C atoms) appears to be an important parameter and can be estimated using the relation proposed by May et al [START_REF] May | From ultrananocrystalline diamond to single crystal diamond growth in hot filament and microwave plasma-enhanced CVD reactors: a unified model for growth rates and grain sizes[END_REF]:

R = 1 1 + 0.3 exp(3430/T S ) + 0.1 exp(-4420/T s )[H 2 ]/[H] (1)
where T s is the surface temperature and [H 2 ] and [H] the gas-phase densities (in cm -3 ) of molecular and atomic hydrogen in the vicinity of the surface, respectively.

Operating conditions for growing diamond on a (100) face, i.e. a substrate temperature of around T s = 1120 K, a methane percentage in H 2 -CH 4 mixture of 5 % to 7 %, a power higher than 3 kW and a pressure from 200 to 400 hPa [START_REF] Achard | The control of growth parameters in the synthesis of high-quality single crystalline diamond by CVD[END_REF][START_REF] Achard | High quality MPACVD diamond single crystal growth: high microwave power density regime[END_REF][START_REF] Gicquel | Quantitative analysis of diamond deposition reactor efficiency[END_REF] give rise to typical densities of [H 2 ] = 10 18 cm -3 and

[H] = 5 • 10 15 cm -3 at the surface [START_REF] Gicquel | Quantitative analysis of diamond deposition reactor efficiency[END_REF], it comes a fraction R of 12.5% at T s = 1120 K, 15% at T s = 1200 K and 20% at T s = 1500 K.

In this article, we will assume equilibrium between the two processes (a)

and (b) and then a constant fraction of radical site at the surface during simulation in order to analyze the interaction of CH 3 radicals with the diamond surface (c).

While molecular dynamics (MD) simulations are well suited for examining basic mechanisms of diamond growth, a rather few ones have been dedicated

to gas phase radicals sticking [START_REF] Alfonso | Molecular-dynamics simulations of methyl-radical deposition on diamond (100) surfaces[END_REF][START_REF] Alman | Molecular dynamics simulation of hydrocarbon reflection and dissociation coefficients from fusion-relevant carbon surfaces[END_REF][START_REF] Eckert | Molecular dynamics simulations of the sticking and etch behavior of various growth species of (Ultra)Nanocrystalline diamond films[END_REF][START_REF] Eckert | On the reaction behaviour of hydrocarbon species at diamond (100) and (111) surfaces: a molecular dynamics investigation[END_REF][START_REF] Huang | Energy dependence of methyl-radical adsorption on diamond (001)-(21) surface[END_REF][START_REF] Kaukonen | Molecular-dynamics simulation of the growth of diamondlike films by energetic carbon-atom beams[END_REF][START_REF] Neyts | Molecular dynamics simulations for the growth of diamondlike carbon films from low kinetic energy species[END_REF][START_REF] Neyts | Influence of internal energy and impact angle on the sticking behaviour of reactive radicals in thin a-c:h film growth: a molecular dynamics study[END_REF][START_REF] Keudell | Surface loss probabilities of hydrocarbon radicals on amorphous hydrogenated carbon film surfaces: Consequences for the formation of re-deposited layers in fusion experiments[END_REF][START_REF] Tichmann | Determination of the sticking coefficient of energetic hydrocarbon molecules by molecular dynamics[END_REF][START_REF] Traskelin | Molecular dynamics simulations of CH3 sticking on carbon surfaces[END_REF][START_REF] Trskelin | Molecular dynamics simulations of CH3 sticking on carbon surfaces, angular and energy dependence[END_REF].

Among them, Eckert et al [START_REF] Eckert | Molecular dynamics simulations of the sticking and etch behavior of various growth species of (Ultra)Nanocrystalline diamond films[END_REF][START_REF] Eckert | On the reaction behaviour of hydrocarbon species at diamond (100) and (111) surfaces: a molecular dynamics investigation[END_REF] have simulated the interaction of various radicals, chosen consistently with experimental conditions and findings of May et al [START_REF] May | Diamond thin films: a 21st-century material[END_REF]: incoming radicals are equilibrated at 2120 K prior releasing to the surface, and simulations are performed for substrate temperatures T s of 800 and 1100 K. CH 3 sticking coefficient is shown to be independent of T s and close to 0.1. However the conditions used in these simulations are relevant for a non-hydrogenated nanodiamond surface and not for hydrogenated single crystalline (100) diamond surface as considered in the present work. A closer work to ours [START_REF] Neyts | Reaction mechanisms of lowkinetic energy hydrocarbon radicals on typical hydrogenated amorphous carbon (a-c:h) sites: A molecular dynamics study[END_REF] is related to the interaction of various radicals, among them on various surface sites, including hydrogenated sites, but for a-C:H

surface not for hydrogenated diamond one. Träskelin et al [START_REF] Traskelin | Molecular dynamics simulations of CH3 sticking on carbon surfaces[END_REF][START_REF] Trskelin | Molecular dynamics simulations of CH3 sticking on carbon surfaces, angular and energy dependence[END_REF] considered interaction with (100) hydrogenated diamond surface and claimed that CH 3 sticking coefficient is less than 0.01. But in this case, the diamond surface is maintained at T s = 0 K, which is far from our conditions. Other studies calculate sticking coefficients ranging from 0.1 to 0.6 [START_REF] Huang | Energy dependence of methyl-radical adsorption on diamond (001)-(21) surface[END_REF][START_REF] Kaukonen | Molecular-dynamics simulation of the growth of diamondlike films by energetic carbon-atom beams[END_REF][START_REF] Tichmann | Determination of the sticking coefficient of energetic hydrocarbon molecules by molecular dynamics[END_REF]. But the calculations were carried out using CH 3 or C atoms at hyperthermal energies, above 1 eV and rising up to 150 eV. This is also far from the conditions explored here, i.e. around 0.1 eV.

The present work is thus intended to determine the CH 3 sticking probabilities on (100) and (111) diamond surfaces in conditions consistent with the microwave plasma assisted growth described in ref. [START_REF] Gicquel | Quantitative analysis of diamond deposition reactor efficiency[END_REF][START_REF] Lombardi | Study of an h 2/CH4 moderate pressure microwave plasma used for diamond deposition: modelling and IR tuneable diode laser diagnostic[END_REF][START_REF] Lombardi | Quantitative detection of methyl radicals in non-equilibrium plasmas: a comparative study[END_REF]. Next section is devoted to the MD procedure consistent with these experimental conditions, while the third section will discuss results for sticking calculation on C(100) and (111) surface at various gas and surface temperature.

MD is a computer simulation technique where the time evolution of a set of interacting atoms is followed by integrating the Newton equations of motion.

The reactive empirical bond order (REBO) potential from Brenner [START_REF] Brenner | Empirical potential for hydrocarbons for use in simulating the chemical vapor deposition of diamond films[END_REF][START_REF] Brenner | A second-generation reactive empirical bond order (REBO) potential energy expression for hydrocarbons[END_REF] which is a semi-empirical potential (analytical function for the interaction model whose fitting parameters are chosen in order to reproduce experimental or ab initio results) is used. The bond order term which appears in the expression of bond energy modulates the bond strength by the local environment (number of bonds, angle...). Therefore chemical reactions can be accounted for, that is to say this potential can handle bond breaking and formation with appropriate change in hybridization. Solving Newton equations of motion requires the knowledge of the initial positions and velocities of all the atoms of the system. The system is split in two subsystems: the diamond substrate and the CH 3 radical. The initial molecular configuration is shown on Fig. 1, where the larger spheres in the lattice represent carbon atoms whereas the smaller spheres at the surface are hydrogen atoms.

Previous studies confirmed that during diamond MW-PACVD growth, the (100) diamond surface is 2×1 reconstructed and partly hydrogenated, depending on the substrate temperature (see Eq. ( 1)). We use C(100) 2×1:H@x, x=100-R standing for x % of hydrogenated 2 × 1 reconstructed (100) diamond surface. The simulated substrate is built using the diamond lattice parameter at the right temperature which is calculated using the relation δa/a = α(T s )δT s where a, T s and α(T s ) are the lattice parameter, the Frauenheim et al [START_REF] Th | Stability, reconstruction, and electronic properties of diamond (100) and (111) surfaces[END_REF] who obtained equilibrium bond length and angles at surface using ab initio calculations and which were confirmed by atomic force microscopy observations [START_REF] Nimmrich | Atomic-resolution imaging of clean and hydrogenterminated c(100)-(2x1) diamond surfaces using noncontact AFM[END_REF].

The surface radical sites (dangling bonds), created by abstraction reaction of surface H atoms by gas-phase impacting atomic H, can be of two main types [START_REF] Skokov | Elementary reaction mechanism for growth of diamond (100) surfaces from methyl radicals[END_REF]: monoradical site (a single dangling bond on a surface carbon)

and biradical sites (two surface radical sites on adjacent carbons). Different variants of these mono-and biradical sites exist depending upon the local surface geometry. From Mankelevich et al [START_REF] Yu | New insights into the mechanism of CVD diamond growth: Single crystal diamond in MW PECVD reactors[END_REF], for typical diamond CVD conditions, the fraction of available biradical sites of all types is ∼ 10 times lower than that of the monoradical sites. We will then assume only monoradical sites at the surface. They are randomly distributed over the surface.

The initial configuration of the CH 3 radical in the gas phase is shown on Massively Parallel Simulator) molecular dynamics solver [START_REF] Plimpton | Fast parallel algorithms for short-range molecular dynamics[END_REF] with a time step of 1 fs.

Equilibration

Before running calculation for the determination of the sticking coefficient of CH 3 , preliminary molecular dynamics simulation are carried out in order to equilibrate both systems (substrate and radical) at the chosen temperature.

The goal of the equilibration phase is not to keep the temperature constant, but rather to ensure that the average temperature of each system be the desired one. The thermalization of the substrate is reached after 5 ps.

The temperature fluctuates around the prescribe experimental value T s = 1120 K. The equilibrated configuration is then stored for the subsequent simulation of the interaction with the incoming CH 3 radical. Concerning the methyl radical, to prevent any bias in the calculation of the velocity distribution, the equilibration is carried out for a set of 8000 CH 3 radicals using a Nosé -Hoover thermostat chain [START_REF] Martyna | Noséhoover chains: The canonical ensemble via continuous dynamics[END_REF]. Once the equilibrium at the targeted gas temperature is reached, the computation continues for successively storing 8000 radicals to make a list of CH 3 configurations which will be used as initial configurations for the subsequent simulation of the interactions with the diamond surface.

Once, the two subsystems are correctly thermalized, they are "merged" in order to calculate the interaction of the methyl radical with the diamond surface. The sticking coefficient describes the probability of an incoming particle to be trapped in a potential well above a surface. Thus, its determination required a statistical study. Before each elementary calculation, an incident radical is picked up from the previously build CH 3 configuration list and is initially placed at a randomly chosen position in the x-and y-directions and at a distance above the substrate located slightly beyond the potential cutoff distance, where interactions cancel each other out. Only selected radicals having a z-velocity component directed toward the surface are released towards the surface. The calculation time for one interaction is inversely related to the magnitude of the z-component of the CH 3 velocity,

v CH 3 z
and thus can be extremely high for the small values of the distribution.

To prevent huge calculation time for sampling the whole velocity distribution, we first assume that, due to steric effect, the contribution of the configurations with v CH 3 z < 100 m/s is negligible. Present calculations show that no methyl radical with vertical velocity component less than 500 m.s -1 sticks (see discussion below). Secondly, the velocity range is split into intervals I v of 200 m.s -1 wide. For each such I v a specific and optimized calculation time is chosen in such a way that each impact is monitored during 1.5 ps.

For each I v interval, 10 000 elementary interactions are calculated. So with 12 intervals, it means that 120000 impacts are simulated for determining the sticking coefficient SC. SC Iv is the sticking coefficient calculated for each interval as:

SC Iv = number of CH 3 with v z ∈ I v which sticks total number of incident CH 3 (2)
Then the global sticking coefficient is obtained by weighting each individual sticking coefficient SC Iv by its corresponding area A Iv under the distribution curve, and summing over all intervals:

SC = Iv SC Iv • A Iv . (3)
3 Results and discussion the surface within the same range but having different vibrational and rotational energies due to equilibration of the methyl radical vapor at three different gas temperatures. In the gas temperature range explored here, the rotational and vibrational energies thus do not seem to contribute significantly on SC values while only a slight dependence on the gas temperature is observed (except a significant effect at 2500 ms -1 ) as previously observed [START_REF] Eckert | On the reaction behaviour of hydrocarbon species at diamond (100) and (111) surfaces: a molecular dynamics investigation[END_REF]. The sticking is clearly increasing when increasing kinetic energy, consistently with all previous studies: the sticking coefficient reaches 0.025 on (111) surface and 0.01 on (100) surface. Moreover, at high incoming velocity, the They have found a sticking coefficient of about 8 • 10 -3 at 1120 K with an estimated accuracy of ± 20%. The measured value is one order of magnitude higher than the calculated one. This can be explained by the polycrystalline nature of the diamond surface used in the experiments compared to the perfect (100) and (111) surface used in calculations. The grain boundaries are expected to provide more active area for sticking. Moreover, they showed surface temperature dependence of recombination efficiency of methyl radicals on diamond surface which confirms our calculations. This is also consistent with another previous work leading to SC < 0.01 [START_REF] Keudell | Surface loss probabilities of hydrocarbon radicals on amorphous hydrogenated carbon film surfaces: Consequences for the formation of re-deposited layers in fusion experiments[END_REF].

Conclusion

MD simulations were carried out for determining the sticking coefficient of CH 3 radical in conditions consistent to MWPACVD growth of thick diamond films. The average sticking coefficient is found to be very low and ranges in [0.28 -0.68] • 10 -3 for the studied conditions of temperatures and substrate orientation, which is consistent with experimental findings and some MD results on hydrogenated diamond surface. Both temperature and velocity range dependence of the sticking coefficient evolution are consistent with previous MD calculations whatever the substrate is.

Table of Contents

Efficient processing of diamond single crystal deposition by moderate pressure and power microwave plasma requires knowledge of important basic data, among them the sticking coefficients of the main deposition precursors. This work provides molecular dynamics simulations for a rigorous of the methyl radical sticking coefficient on hydrognated diamond single crystal substrate.

Figure 1 :

 1 Figure 1: Initial molecular configuration of the diamond substrate

Fig. 2 .

 2 Fig. 2. The initial C-H bond length and H-C-H bond angle are taken from Marynick et al [42] who reproduced, through ab initio calculations, the planar structure of the methyl radical. The initial velocities of all the atoms

Figure 2 :Figure 3 :

 23 Figure 2: Initial molecular configuration of the methyl radical

Fig. 4

 4 Fig. 4 displays the sticking coefficient for each velocity range at the three different gas and substrate temperatures (1120 K, 1200 K and 1500 K) and for both diamond surfaces, C(100) and C(111). The symbols located on a vertical line represent values of sticking coefficient of CH 3 onto both (100) and (111) surfaces, for CH 3 having translational kinetic energy normal to

Figure 4 :Figure 5 :

 45 Figure 4: Sticking coefficients of CH 3 onto diamond surfaces versus the velocity range

http://lammps.sandia.gov

Acknowledgments

ANR is gratefully acknowledged for grant # ANR2011BS09029 PlasBorD-DiAM. ALINEOS SA is gratefully acknowledged for configuring the multiprocessor high performance computer used in this study.