Loïc Peletan
email: loic.peletan@insa-lyon.fr

Sébastien Baguet
email: sebastien.baguet@insa-lyon.fr

Mohamed Torkhani
email: mohamed.torkhani@edf.fr

Georges Jacquet-Richardet

Quasi-Periodic Harmonic Balance Method for Rubbing Self-Induced Vibrations in Rotor-Stator Dynamics

Keywords: quasi-periodic harmonic balance method, self-induced excitation, automatic harmonic selection, nonlinear rotordynamics, rotor-stator contact

A quasi-periodic Harmonic Balance Method (HBM) coupled with a pseudo-arc length continuation algorithm is developed and used for the prediction of the steady state dynamic behaviour of rotor-stator contact problems. Quasiperiodic phenomena generally involve two incommensurable fundamental frequencies and at present the Harmonic Balance Method has been adapted to deal with cases where those frequencies are known. The problem here is to improve the procedure in order to be able to deal with cases where one of the two fundamental frequencies is a priori unknown, in order to be able to reproduce self-excited phenomena such as the so-called quasi-periodic partial rub. Considering the proposed developments, the unknown fundamental frequency is automatically determined during calculation and an automatic harmonic selection procedure gives both accuracy and performance improvements. The application is based on a Jeffcott rotor model and results obtained are compared with traditional time marching solutions. The modified quasi-periodic HBM appears one order of magnitude faster than transient simulations while providing very accurate results.

Introduction

Improving turbomachines efficiency may be achieved by reducing the operating gap between rotating and stationary parts (rotor and stator). However this leads to an increased risk of rotor-to-stator contact and, for safety reasons, designers must ensure that such interactions cannot have serious consequences. When an accurate prediction of the system response is required, it is not sufficient to consider linear models. When nonlinearities are considered, highly complex dynamic behaviours are possible, even in the case of very simple rotors [START_REF] Jacquet-Richardet | Rotor to stator contacts in turbomachines. review and application[END_REF]. The simplest rotor model encountered is the Jeffcott model described by only one node and two degrees of translational freedom (one for each radial direction). Edwards [START_REF] Edwards | The influence of torsion on rotor/stator contact in rotating machinery[END_REF] studied a Jeffcott rotor with an additional torsional degree of freedom, when subjected to a mass unbalance and contact with a circular, rigid, fixed stator. He modelled the contact using a penalty method with a Coulomb friction law and showed that even when the rotor is only excited by a periodic force, this simple system may exhibit a rich variety of phenomena due to the presence of nonlinearity.

Many authors have studied nonlinear rotor models, using time marching techniques [START_REF] Edwards | The influence of torsion on rotor/stator contact in rotating machinery[END_REF][START_REF] Feng | Rubbing phenomena in rotor-stator contact[END_REF][START_REF] Popprath | Nonlinear dynamics of a rotor contacting an elastically suspended stator[END_REF][START_REF] Han | Periodic motions of a dual-disc rotor system with rub-impact at fixed limiter[END_REF][START_REF] An | Dynamic response of a rub-impact rotor system under axial thrust[END_REF]. This approach, although essential when investigating the system's transient behaviour [START_REF] Feng | Rubbing phenomena in rotor-stator contact[END_REF][START_REF] Grapis | Overcritical high-speed rotor systems, full annular rub and accident[END_REF][START_REF] Braut | Dynamic analysis of the rotor-stator contact due to blade loss[END_REF], is not optimal in terms of CPU time when the steady state response is investigated, because many cycles need to be simulated before the transient behaviour disappears completely. This can be problematic when parametric studies need to be carried out and/or when more complex models are used (such as finite element rotor models [START_REF] Han | Periodic motions of a dual-disc rotor system with rub-impact at fixed limiter[END_REF][START_REF] Qin | Nonlinear responses of a rub-impact overhung rotor[END_REF]).

On the other hand, some authors such as Jiang [START_REF] Jiang | Stability analysis of sliding whirl in a nonlinear jeffcott rotor with cross-coupling stiffness coefficients[END_REF][START_REF] Jiang | Determination of the global responses characteristics of a piecewise smooth dynamical system with contact[END_REF] or Bently et al. [START_REF] Bently | Full annular rub in mechanical seals, part ii: Analytical study[END_REF] use analytical tools to study simple rotors and determine the limits of existence for different types of behaviours, according to a set of dimensionless parameters. Considering Jiang's model the stator is circular, rigid and Fig. 1 Schematic representation of a Jeffcott rotor with a circular rigid stator, after [START_REF] Jiang | Determination of the global responses characteristics of a piecewise smooth dynamical system with contact[END_REF].

fixed and contact is modelled by a penalty law using a dry friction Coulomb approach (Fig. 1). His model is able to exhibit the four classical types of steady-state rotor behaviour, resulting from unbalance and rotor-to-stator contact. These four types of behaviour are:

-'no-rub motion' in which there is no contact between the rotor and the stator. The orbit is circular and the motion is periodic. -'full annular rub' in which the rotor remains in permanent contact with the stator. The orbit is also circular and the motion periodic. -'partial rub' in which the rotor intermittently comes into contact with the stator. The motion is quasi-periodic. -'backward whip' in which the rotor rolls (with slippage) on the inner surface of the stator. This motion is selfexcited and its frequency is not necessarily a direct function of rotational speed. The phenomenon is violent and a priori quasi-periodic.

The governing equations of the rotor/stator contact system can be written as

X + 2ξX + βX + Θ Å 1 - R 0 R ã X -µsign (V rel ) Y = Ω 2 cos Ωτ Y + 2ξY + βY + Θ Å 1 - R 0 R ã µsign (V rel ) X + Y (1)
= Ω 2 sin Ωτ

V rel = R disk Ω + Rω b
Jiang's study [START_REF] Jiang | Determination of the global responses characteristics of a piecewise smooth dynamical system with contact[END_REF] is based on the following dimensionless variables: 2 Rotor response characteristics in the Ω-µ plane with β = 0.04, ξ = 0.05, R 0 = 1.05 and R disk = 20R 0 , after [START_REF] Jiang | Determination of the global responses characteristics of a piecewise smooth dynamical system with contact[END_REF]. Some characteristic values of Ω associated to µ = 0.2 are highlighted.

X = x e , Y = y e , R 0 = r 0 e , R disk = r disk e , V rel = v rel e , 2ξ = c √ k b m , β = k s k b , ω s = » k b m , Ω = ω ω s , ω b = ω w ω s , Ω DF =0.
Ω u Ω l Fig.
where x and y represent the radial displacements along the x and y axes respectively, e is the mass eccentricity, r 0 the initial clearance between the rotor and the stator, r disk the disk radius, v rel the relative velocity between the rotor and the stator, c the damping factor, m the mass of the rotor, k s and k b are the stiffnesses of the rotor and stator respectively, ω is the rotational speed of the rotor, ω w the whirl frequency of the rotor and ω s the natural frequency of the rotor system without clearance. The coefficient of friction is denoted µ.

Fig. 2 is extracted from Jiang's paper [START_REF] Jiang | Determination of the global responses characteristics of a piecewise smooth dynamical system with contact[END_REF]. It shows the analytical prediction of the response characteristics of the rotor in the Ω-µ plane. The straight lines Ω l and Ω u indicate the rotational speed limits where rotor eccentricity reaches the initial clearance. Thus, the no-rub motion can appear only at rotational speeds below Ω l and above Ω u . The curves HP and SN were obtained numerically and correspond to the limits of full annular rub motion. To the right of line HP, the full annular rub becomes unstable and motion becomes quasi-periodic (Hopf bifurcation). To the right of line SN, the full annular rub no longer exists (saddle node bifurcation). Line DF is the backward whip limit and line DW indicates the point beyond which backward whip is triggered by unbalance. Some of Jiang's analytical equations had to be solved with the help of numerical calculations. It explains the dots on the HP and SN curves.

There are some regions in which several behaviours coexist. Within regions labelled '0' no-rub and dry whip may exist. In the region labelled '1', both full annular rub and dry friction backward whirl are possible. The region labelled '2' shows the presence of partial rub and backward whirl. In '3', full annular rub and no-rub motions are simultaneously present. Finally, a small region '4' exhibits no-rub, full annular rub and backward whip possibilities. Time integration methods would lead to only one of these solutions, depending on initial conditions. Analytical methods have been extended to multi-mode rotor models by Childs and Bhattacharya [START_REF] Childs | Prediction of dry-friction whirl and whip between a rotor and a stator[END_REF]. Their predictions appear to be in good agreement with numerical simulations and experiments. Nevertheless, these analytical methods require non-trivial mathematical developments and their applicability is limited when considering complex models.

Alternatively, although widely used to solve other types of problems, numerical techniques such as shooting or the Harmonic Balance Methods (HBM) are rarely used for rotorstator contact problems. Although such techniques do not require long transient simulations and are in practice much faster, they are usually limited to the study of periodic motions. This limitation explains that those techniques are generally restricted to the study of geometric nonlinearities or oil-film bearing nonlinearities [START_REF] Sundararajan | An algorithm for response and stability of large order non-linear systems -application to rotor systems[END_REF], because with this kind of nonlinearities the system response remains mostly periodic. Fortunately, extension of the HBM to quasi-periodicity is possible by introducing the concept of multidimensional time and has already been performed for various applications. Lau [START_REF] Lau | Incremental harmonic balance method with multiple time scales for aperiodic vibration of nonlinear systems[END_REF] studied the quasi-periodic free response of a clamped beam with large displacements. Legrand [START_REF] Legrand | Modèles de prédiction de l'intéraction rotor/stator dans un moteur d'avion[END_REF][START_REF] Legrand | ndimensional harmonic balance method extended to non-explicit nonlinearities[END_REF] used the quasi-periodic HBM for the computation of limit cycles for autonomous systems with application to aircraft turbomachinery with rotor-stator interaction. Coudeyras [START_REF] Coudeyras | Periodic and quasi-periodic solutions for multi-instabilities involved in brake squeal[END_REF] used a similar approach applied to the study of break squeal. Complex nonlinear normal modes can be computed with the help of quasi-periodic HBM as done by Laxalde and Thouverez [START_REF] Laxalde | Complex non-linear modal analysis for mechanical systems: Application to turbomachinery bladings with friction interfaces[END_REF] for blade-casing contact problems. The quasiperiodic HBM coupled with an arc-length continuation scheme was also used by Guskov [START_REF] Guskov | Multi-dimensional harmonic balance applied to rotor dynamics[END_REF] for the forced response of a multiple shaft rotor subject to quasi-periodic multi-excitations with known frequencies. Some variants can be found in the literature. For example, instead of computing a quasi-periodic trajectory which takes place on the surface of an invariant torus, Schilder et al. [START_REF] Schilder | Fourier methods for quasi-periodic oscillations[END_REF] focus on the direct computation of the torus itself and solve a so-called Invariance PDE that a torus function must satisfy, with application to non-linear electrical engineering.

In the following, the nonlinear Jeffcott rotor model described in [START_REF] Jiang | Determination of the global responses characteristics of a piecewise smooth dynamical system with contact[END_REF] is first considered and analysed using a classical time marching procedure. The goal is to simulate the transient rotor behaviour in order to analyse carefully the responses obtained. The highlighted important results lead to the proposition of an improved Harmonic Balance Method associated with a pseudo arc-length continuation technique, able to reproduce efficiently quasi-periodic phenomena where nothing is a priori known concerning the frequency induced by the contact mechanism.

Nonlinear Jeffcott rotor -numerical approach

The results shown in Fig. 3 have been examined via numerical simulations based on a time marching procedure.

Considering the same parameter values as those presented in Fig. 2, simulations were run with a fixed value for the friction coefficient µ and varying values of Ω. Solutions were calculated for discrete values of Ω and the initial conditions used at each step were given by the final steady state solution of the previous one (incremental continuation technique). For each value of Ω, 100 revolutions were simulated in order to reach the steady-state behaviour. Twenty-five additional revolutions were then simulated and saved for analysis. This procedure was performed twice for each value of µ. In the first round, simulations were made with progressively greater values of Ω, whereas in the second progressively smaller values of Ω were used. As numerical simulations require the use of non dimensionless variables, the following parameters were set arbitrarily: e = 0.1 m, k s = 100 N.m -1 , m = 1 kg and the other parameters deduced from those. All the simulations were performed using Code Aster [START_REF] Edf R&d | Code aster: A general code for structural dynamics simulation under gnu gpl licence[END_REF] and user scripts in Scipy [START_REF] Jones | Scipy: Open source scientific tools for Python[END_REF] environments. Although simulations involved many values of µ, only the results obtained for µ = 0.2, which are representative of all possible solutions, are presented in the following. The response curve obtained for µ = 0.2, β = 0.04, ξ = 0.05, R 0 = 1.05 and R disk = 20R 0 is shown in Fig. 3. The y-axis represents rotor eccentricity normalized with the initial clearance, on a logarithmic scale, and the x-axis is the normalized rotational speed Ω = ω ω s . Different behaviours are obtained according to the analytical prediction (Fig. 2).

During the sweep-up, the no-rub motion is first observed (Fig. 4). At Ω l ≈ 0.16, the full annular rub regime starts with a rotor eccentricity slightly larger than the initial gap due to the contact elasticity (Fig. 5). The transition between no-rub and full annular rub (which corresponds to the Ω l line in Fig. 2) occurs between approximately Ω = 0.155 and Ω = 0.160, which is highly consistent with the analytical value (Ω l = 0.1536). Beyond Ω HP ≈ 0.29, rotor motion becomes quasi-periodic. The orbit at Ω = 0.33 is shown in Fig. 6 and illustrates the partial rub motion predicted, where the rotor does not remain in permanent contact with the stator and the maximum rotor eccentricity is greater than that predicted during full annular rub. At Ω DW ≈ 0.46, eccentricity suddenly increases by a factor of ten. Although slightly different from the value found analytically (≈ 0.41), the value obtained is coherent. In this backward dry whip regime, rotor's orbit is almost circular, as shown in Fig. 7. The rotor then remains in the backward whip regime until the end of the sweep-up.

During the sweep-down, the rotor starts with the no-rub regime, until Ω u ≈ 0.86 where a sudden increase in eccentricity occurs, transferring the motion of the rotor to the backward whip regime. The rotor remains in this regime until a very small value of Ω DF ≈ 0.04 is reached, at which eccentricity suddenly drops leading to the no rub regime. This value of Ω DF is also in good agreement with the analytical value (Fig. 2). Transient simulations are used to analyse more deeply the nature of the quasi-periodic regimes. A Fast Fourier Transform of the nonlinear contact force, during partial rub at Ω = 0.33, is shown Fig.

Ω 0a =1 Ω 0b =2.08 Ω 1a =4.08 Ω 1b =5.16 Ω 2a =7.16 Ω 2b =8.24 Ω 3a =10.2 Ω 3b =11.
3 those results, it appears that the frequencies obtained are all linked by a linear combination of two incommensurable frequencies. The lowest frequency Ω 1 is associated to rotation (1 × Ω) while the second Ω 2 cannot be determined a priori since it is not commensurable with Ω 1 , here Ω 2 ≈ 2.08 × Ω.

Ω na =(n+1)Ω 1 + (n+0)Ω 2 Ω nb =(n+0)Ω 1 + (n+1)Ω 2 with Ω 1 =1
For every simulated partial rub motions, spectral peaks appear in the contact force power spectrum at linear combinations of the two incommensurable frequencies, according to the following relations:

Ω na = (n + 1)Ω 1 + (n + 0)Ω 2 for n ∈ N (2)
Ω nb = (n + 0)Ω 1 + (n + 1)Ω 2 for n ∈ N (3)
Peaks appear in pairs and decrease in amplitude with n (Fig. 8). Values of Ω na and Ω nb are given in Table 1 and are coherent with those highlighted in Fig. 8.

Spectral analysis of the contact force during partial rub for several values of Ω shows not only that Ω 1 remains equal to the rotating frequency (Ω 1 = 1 × Ω), but also that Ω 2 is related to Ω 1 without any trivial relationship (Ω 2 = f (Ω 1 )). The only evident point is that Ω 2 decreases when Ω increases. The evolution of Ω 2 with Ω 1 is given in Fig. 9. A nonlinear decrease of Ω 2 with Ω 1 is observed. As the two fundamental frequencies are incommensurable and due to limited spectral resolution, only approximate measurements of Ω 2 can be obtained, with lower and upper bounds indicated by error bars in Fig. 9.

Those direct simulations are highly computer time demanding. On the other hand the harmonic balance method is very efficient but is limited to the study of periodic solutions or quasi-periodic solutions in the case where the relationship between the different frequency components is known. In the following, after a brief presentation of the periodic HBM and the pseudo arc-length continuation method, an extension of the HBM is proposed. As required here, the method handles situations where one of the two fundamental frequencies is a priori unknown.

3 Periodic Harmonic Balance Method (HBM)

Basics of the HBM

The classical Harmonic Balance Method (HBM) is used to compute periodic solutions of dynamical systems. With this method the equations of motion are solved in the frequency domain, rather than in the time domain. The general equation for dynamical systems is:

M q(t) + C q(t) + Kq(t) + f(q(t), q(t)) -p(t) = 0 ( 4
)
where q is the displacement vector for all of the n degrees of freedom; K, C and M are the generalized n × n stiffness, damping and mass matrices; f is the nonlinear force vector and p the external excitation force vector.

When the external excitations are periodic, it is possible to assume that a steady state solution for Eq. ( 4) exists, and that this solution is also periodic. The displacements, the external and the nonlinear forces can thus be written as truncated Fourier series with N harmonics:

q(t) = Q 0 + N j=1 Q 2 j-1 cos( jωt) + Q 2 j sin( jωt)
Not all the harmonics contribute to the solution in the same way. Some harmonics are prevailing while some others have a nil contribution. If only the harmonics of interest (see section 4.3 for harmonic selection) are retained in the solution, then the truncated Fourier series for displacements, external and nonlinear forces become:

q(t) = Q 0 + N j=1 Q 2 j-1 cos(m j ωt) + Q 2 j sin(m j ωt) (5)
p(t) = P 0 + N j=1 P 2 j-1 cos(m j ωt) + P 2 j sin(m j ωt) (6)
f(t) = F 0 + N j=1 F 2 j-1 cos(m j ωt) + F 2 j sin(m j ωt) (7)
Q = Q 0 T , Q 1 T , . . . , Q T 2N T , F = F T 0 , F T 1 , . . . , F T

2N

T and P = P T 0 , P T 1 , . . . , P T

2N

T are vectors of the Fourier coefficients for displacements, nonlinear forces and external excitations, respectively. ω is the fundamental frequency of the external excitation, and m j ∈ N (1 ≤ j ≤ N) represent the indexes of the N harmonics of interest retained in the solution.

As described in [START_REF] Petrov | Analytical formulation of friction interface elements for analysis of nonlinear multi-harmonic vibrations of bladed disks[END_REF] for example, Eqs. (5) ( 6) and ( 7) can be substituted into Eq. ( 4) and a Galerkin procedure applied to transform the nonlinear differential Eq. (4) of dimension n, into an algebraic nonlinear system of dimension n HBM = n(2N + 1):

R(Q, ω) = Z(ω)Q + F(Q) -P = 0 (8)
where Z(ω) = diag K, Z 1 (ω), . . . , Z j (ω), . . . , Z N (ω) and:

Z j (ω) = ñ K -m j 2 ω 2 M m j ωC -m j ωC K -m j 2 ω 2 M ô
Eq. ( 8) has to be solved for Q. As this equation is still nonlinear, an incremental-iterative Newton-Raphson method is used to derive a correct solution.

The Alternating Frequency Time (AFT) algorithm [START_REF] Cameron | An alternating frequency/time domain method for calculating the steady-state response of nonlinear dynamic systems[END_REF] gives an efficient way to compute the nonlinear term F(Q) and its derivative ∂F(Q)/∂Q at each Newton-Raphson iteration. The nonlinear forces are usually much easily evaluated in the time domain than in the frequency domain. Moreover, the AFT scheme uses fast direct and inverse Fourier transforms to quickly compute the nonlinear forces in the time domain and then switch back to the frequency domain. A combination with a pseudo-arc length continuation [START_REF] Sundararajan | An algorithm for response and stability of large order non-linear systems -application to rotor systems[END_REF] allows following efficiently the solution branches. The stability of the solutions is determined using the Floquet's theory [START_REF] Peletan | A comparison of stability computational methods for periodic solution of nonlinear problems with application to rotordynamics[END_REF].

Pseudo arc-length continuation method

The interest of the pseudo arc-length continuation method lies in its ability to follow response branches beyond turning points and thus completely describe a complex solution branch with potentially several solutions for a given value of ω. Let M (i) Q (i) , ω (i) be a point that satisfies R(Q (i) , ω (i) ) = 0. Then the goal is to obtain a point M (i+1) Q (i+1) , ω (i+1) that also satisfies R(Q (i+1) , ω (i+1) ) = 0, with Q (i+1) = Q (i) + ∆Q (i) and ω (i+1) = ω (i) + ∆ω (i) . The continuation procedure consists in two steps: a prediction and corrections.

ω Q e (i) e (i+1)
Q (i) ,ω (i) Q (i+1) ,ω (i+1) (0) Q (i+1) , (0) ω (i+1) (k) Q (i+1) , (k) ω (i+1)

Prediction

During the prediction step, the solution for the next value of ω is assessed. The prediction is made in a direction tangent to the solution branch. The following first order approximation can be made:

R(Q (i+1) , ω (i+1) ) =0 ≈ R(Q (i) , ω (i) ) =0 +∆Q (i) ∂R ∂Q M (i) + ∆ω (i) ∂R ∂ω M (i) (9)
with ∂R ∂Q = Z(ω) + dF dQ and ∂R ∂ω = dZ dω Q. The norm of the tangent vector e (i) = [∆Q T (i) , ∆ω (i) ] is:

e (i) 2 = ∆Q T (i) ∆Q (i) + ∆ω 2 (i) Let consider a (i) = ∆ω (i) , ∆Q (i) = a (i) ∆ Q(i)
and the norm of the tangent vector set at unity, then:

a (i) = ± 1 » ∆ QT (i) ∆ Q(i) + 1 (10)
The sign of a (i) is chosen so that two consecutive tangent vectors have a positive scalar product (direction conservation). Eq. ( 9) becomes:

∆Q (i) ∂R ∂Q M (i) + ∆ω (i) ∂R ∂ω M (i) = 0 (11)
Considering the expressions of a (i) (a (i) = ∆ω (i) ) and ∆Q (i) (∆Q (i) = a (i) ∆ Q(i) ), Eq. (11) becomes:

ä(i) ∆ Q(i) ∂R ∂Q M (i) + ä(i) ∂R ∂ω M (i) = 0 (12)
Hence, the procedure to compute e (i) can be summarized as follows:

-

Solve ∂R ∂Q M (i) ∆ Q(i) = -∂R ∂ω M (i) for ∆ Q(i) -Compute a (i) from Eq. (10) -Deduce e (i) = [∆Q T (i) , a (i) ]
Once the normalised tangent vector is calculated, prediction can be made:

Q (0) (i+1) = Q (i) + ∆s (i) • e (i) (13)
with ∆s (i) being the step size.

Produit par HAL -16 Sep 2014

Corrections

The prediction does not a priori satisfies equilibrium. It has to be corrected iteratively in order to cancel the error. The k th iteration gives:

R(Q (k) (i+1) , ω (k) (i+1) ) =R (k) = Z(ω (k) (i+1) ) =Z (k) Q (k) (i+1) + F(Q (k) (i+1) ) =F (k)
-P

The (i + 1) notation will now be dropped for clarity. Corrections are forced to be orthogonal to the tangent vector e. Thus, the following system has to be solved: ñ

∂R (k) ∂Q ∂R (k) ∂ω ∆ QT (i) ∆ω (i) ô ï ∆Q (k+1) ∆ω (k+1) ò = ï -R (k) 0 ò ( 14
)
Both ω and Q are corrected at each iteration. Iterations are carried out until convergence. This continuation procedure is illustrated in Fig. 10.

Application to the Jeffcott rotor

The results provided by the HBM with pseudo arc-length continuation and stability analysis have been added to the response curve in Fig. 11. As expected, only no-rub and full annular rub motions are obtained by the periodic HBM procedure. The stability analysis shows a loss of stability of the full annular rub beyond Ω HP ≈ 0.29. The Floquet multipliers of the unstable solution at Ω ≈ 0.3 can be seen in Fig. 12. A pair of Floquet multipliers are not within the unit circle, with a non-zero imaginary component indicating a Hopf bifurcation and a loss of stability. This Hopf bifurcation corresponds to the appearance of a second frequency in the real response of the rotor, and the initiation of a quasi-periodic regime. This means that although, theoretically, the periodic full annular rub response still exists beyond this point, it is in practice no longer possible to achieve this response, due to its instability. The value of Ω HP obtained with the HBM is very consistent with the predicted value (see Fig. 2).

The HBM is able to predict only a portion of the time transient simulations. However, stability analysis can detect unstable solutions and then the DF curve can be determined in addition to the Ω l and Ω u lines. The DW and DF lines cannot be predicted by the periodic HBM. Nonetheless, a perfect match is found between stable HBM solutions and time integration results. In the following it will be shown that the HBM can be extended to quasi-periodic analysis allowing the partial rub and backward whip regimes to be predicted too.

Quasi-periodic HBM

As shown above, only the external excitation frequency ω 1 is known here and the second fundamental frequency ω 2 remains an additionnal unknown.

Extension of the HBM to bi-periodic solutions

The concept of hyper-time is used to extend the Harmonic Balance Method to quasi-periodicity. When two fundamental frequencies are taken into account, the hyper-time is bidimensional and the solution is supposed to be periodic for each of the dimensions. Thus, the form of the truncated Fourier series presented in Eqs. ( 5), ( 6) and ( 7) must be adapted. Displacement q is approximated by trigonometric polynomials with two basic frequencies ω 1 and ω 2 :

q(t) = Q 0 + N j=1 Q 2 j-1 cos((m j1 ω 1 + m j2 ω 2 )t) +Q 2 j sin((m j1 ω 1 + m j2 ω 2 )t) (15)
where m j1 and m j2 ∈ Z (1 ≤ j ≤ N) represent the indexes of the N harmonics of interest retained in the solution for ω 1 and ω 2 respectively. Introducing the bi-dimensional vector of basic frequencies ω = [ω 1 , ω 2 ], the bi-dimensional harmonic index m j = m j1 , m j2 , and the hyper-time θ = [θ 1 , θ 2 ] = ωt, this equation can be rewritten in a form very Produit par HAL -16 Sep 2014

Q ω 1 ω 2 ω 2 = f(ω 1 ) R(Q,ω 1 ,ω 2 )=0 e 1 e 2 (k)
Δ s similar to Eq. ( 5) used for periodic HBM:

Q (i) ,ω 1(i) ,ω 2(i) Q (i+1) ,ω 1(i+1) ,ω 2
q(θ) = Q 0 + N j=1 Q 2 j-1 cos m j , θ + Q 2 j sin m j , θ (16)
where •, • stands for the classical dot product. The Galerkin procedure in this case leads to the following algebraic system:

R(Q, ω 1 , ω 2 ) = Z(ω 1 , ω 2 )Q + F(Q) -P = 0 ( 17
)
where Z(ω 1 , ω 2 ) = diag K, Z 1 , . . . , Z j , . . . , Z N and:

Z j = ï K -m j , θ 2 M m j , θ C -m j , θ C K -m j , θ 2 M ò
As for periodic HBM, the purpose is to minimise the residual R. The same tools (Newton-Raphson solver, AFT, pseudo arc-length continuation) are used but need to be adapted to the quasi-periodic case. In particular, the pseudo arc-length continuation scheme must be adapted so that the second fundamental frequency is solved during the Newton-Raphson procedure.

Extension of the pseudo arc-length continuation procedure

The pseudo arc-length continuation algorithm is extended to a parametric study with two parameters (ω 1 and ω 2 ). These two parameters are not independant but are linked by an a priori unknown function (see Fig. 13).

As before

M (i) (Q (i) , ω 1(i) , ω 2(i)
) is a point that satisfies R(Q (i) , ω 1(i) , ω 2(i) ) = 0, and we are looking for a point

M (i+1) (Q (i+1) = Q (i) +∆Q (i) , ω 1(i+1) = ω 1(i) +∆ω 1(i) , ω 2(i+1) = ω 2(i) + ∆ω 2(i) ) that also satifies R(Q (i+1) , ω 1(i+1) , ω 2(i+1) ) = 0.

Prediction

Here also, the following first order approximation can be made:

R(Q (i+1) , ω 1(i+1) , ω 2(i+1) ) =0 ≈ R(Q (i) , ω 1(i) , ω 2(i) ) =0 + ∆Q (i) ∂R ∂Q M (i) + ∆ω 1(i) ∂R ∂ω 1 M (i) + ∆ω 2(i) ∂R ∂ω 2 M (i) (18)
with,

∂R ∂Q = Z(ω 1 , ω 2 ) + dF dQ ∂R ∂ω 1 = ∂Z ∂ω 1 Q ∂R ∂ω 2 = ∂Z ∂ω 2 Q The norm of the tangent vector e 1(i) = [∆Q T (i) , ∆ω 1(i) , ∆ω 2(i) ] is: ||e 1(i) || 2 = ∆Q T (i) ∆Q (i) + ∆ω 2 1(i) + ∆ω 2 2(i)
Let's denote c (i) the local relationship between variations of ω 1 and ω 2 so that c (i) = ∆ω 2 ∆ω 1 . As this relationship is unknown, the following finite differences approximation can be made, using the previouly calculated points:

c (i) = ω 2(i) -ω 2(i-1) ω 1(i) -ω 1(i-1)
The developments are totally similar to those previously presented with a (i) = ∆ω 1(i) » 1 + c 2 (i) and ∆Q (i) = a (i) ∆ Q(i) . In this case, the procedure can be summarised as follows:

-Solve the following system for ∆ Q(i) :

∂R ∂Q M (i) ∆ Q(i) = - 1 » 1 + c 2 (i) Ç ∂R ∂ω 1 M (i) + c (i) ∂R ∂ω 2 M (i) å
-Compute a (i) from Eq. ( 10) -Deduce ∆ω 1(i) and ∆ω 2(i) :

∆ω 1(i) = a (i) » 1 + c 2 (i) ; ∆ω 2(i) = c (i) ∆ω 1(i)
This gives vector

e 1 = [∆ Q(i) , ∆ω 1(i) , ∆ω 2(i) ]. The prediction of point M (i+1
) is made at a distance ∆s from M i along the direction given by e 1 :

M (0) (i+1) = M (i) + ∆s.e 1

Corrections

The residual vector R at iteration k is given by:

R(Q (k) (i+1) , ω (k) 1(i+1) , ω (k) 2(i+1) ) =R (k) = Z(ω (k) 1(i+1) , ω (k) 2(i+1) ) =Z (k) Q (k) (i+1) + F(Q (k) (i+1) ) =F (k)
-P

Produit par HAL -16 Sep 2014

As before, subscript (i + 1) will be dropped for clarity. Corrections are made orthogonally to the tangent vector e 1 and the system at iteration k becomes:

ñ ∂R (k) ∂Q ∂R (k) ∂ω 1 ∂R (k) ∂ω 2 ∆ QT (i) ∆ω 1(i) ∆ω 2(i) ô   ∆Q (k+1) ∆ω (k+1) 1 ∆ω (k+1) 2   = ï -R (k) 0 ò ( 19
)
This system is underdetermined because of the presence of the additional unknown ∆ω (k+1)

2

. In fact, since θ 2 = ω 2 t does not appear in the forcing term of the equations of motion (4), if q(θ 1 , θ 2 ) is a solution then a phase-shifted function q(θ 1 , θ 2 + σ) is also a solution, for any value of σ, i.e. the phase is arbitrary. To define a unique solution an extra constraint is needed. This additional equation is often referred to as a "phase condition" and consists in fixing one of the variables other than ω 1 and ω 2 , i.e. one component of vector Q or, within the Newton-Raphson procedure, one component of correction ∆Q. There is no systematic way for choosing the variable to be fixed. Several phase conditions can be found in textbooks dealing with periodic solutions [START_REF] Nayfeh | Applied Nonlinear Dynamics: Analytical, Computational, and Experimental Methods[END_REF][START_REF] Doedel | Lecture notes on numerical analysis of nonlinear equations[END_REF]. They can be considered as "orthogonality conditions" because they can be written as a null scalar product between the vector of unknowns and a vector e 2 with non-null components along ∆Q only, i.e. e 2 = e 2Q , 0, 0 . With this additional equation, the system to solve becomes:

  ∂R (k) ∂Q ∂R (k) ∂ω 1 ∂R (k) ∂ω 2 ∆ QT (i) ∆ω 1(i) ∆ω 2(i) e 2Q 0 0     ∆Q (k+1) ∆ω (k+1) 1 ∆ω (k+1) 2   =   -R (k) 0 0   (20)
Two constructions for e 2 are proposed in the present paper, which differ in the way ∆Q is constrained. The first choice is the simplest one and consists in directly fixing one arbitrarily chosen component of ∆Q by considering the following phase equation:

∂∆q i ∂θ 2 (θ) = 0 (21)
with ∆q i being the displacement correction of the i th degree of freedom with i arbitrarily chosen between 1 and n. This means that a zero phase shift relatively to the second hypertime dimension θ 2 is forced for this particular degree of freedom. Once processed to the frequency domain, Eq. ( 21) becomes:

N j=1 m j2 ∆Q 2 ji = 0 (22)
Using Eq. ( 22) is equivalent to using system [START_REF] Guskov | Multi-dimensional harmonic balance applied to rotor dynamics[END_REF] with the following e 2 vector:

with m j2 at (n(2 j -1) + i) th position, j ∈ [1..N]. With this phase equation, the system is well determined and the Newton-Raphson procedure can be carried out. This specific e 2 vector is used for each iteration of the Newton-Raphson solver. Alternatively, a second possible choice for e 2 relies on geometric considerations. It consists in requiring the current correction to be orthogonal to the previous one as in the conjugate gradient procedure. Thus, vector e 2 can be chosen as e (k) 2 = [∆Q (k)T , 0, 0]. For the applications considered in this paper, this approach appears slightly more effective than the previous one. However, the drawback here is that it is not possible to build vector e 2 at the first iteration of the Newton-Raphson procedure. One way to overcome this problem is to perform a least square minimisation step on the under determined system [START_REF] Laxalde | Complex non-linear modal analysis for mechanical systems: Application to turbomachinery bladings with friction interfaces[END_REF] or to use the phase equation for the first iteration.

Automatic harmonic selection

Usually, simulations are carried out using an arbitrary set of harmonics. For instance, Chua and Ushida in [START_REF] Chua | Algorithms for computing almost periodic steady state response of non-linear systems to multiple input frequencies[END_REF] and Kim and Choi in [START_REF] Kim | A multiple harmonic balance method for the internal resonant vibration of a non-linear jeffcott rotor[END_REF] select harmonics as a function of an arbitrary number M. Harmonics are then chosen so that:

|m j1 | + |m j2 | ≤ M
Legrand in [START_REF] Legrand | Modèles de prédiction de l'intéraction rotor/stator dans un moteur d'avion[END_REF] uses an even less restrictive formula:

|m j1 | ≤ M and |m j2 | ≤ M
Using one of these formulae leads to use a large number of harmonics because the total number of harmonics rapidly increases with M and the number of harmonics has a direct impact on performance. However many of these harmonics may be unnecessary and some important harmonics may not have been selected.

To overcome these problems, a technique for automatic selection of predominant harmonics can be implemented. Several selection criteria were proposed in the context of conventional periodic HBM. Among those, a simple selection procedure based on the error on the approximated nonlinear forces was presented by Laxalde in [START_REF] Laxalde | Etude d'amortisseurs non-linéaires appliqués aux roues aubagées et aux systèmes multi-étages[END_REF]. Another criterion was presented by Jaumouillé in [START_REF] Jaumouille | An adaptive harmonic balance method for predicting the nonlinear dynamic responses of mechanical systemsapplication to bolted structures[END_REF] in order to take into account the global system behavior. It relies on the evolution of strain energy with an increasing number of harmonics. However, this incremental criterion requires an extra calculation cost for strain energy and the selection procedure can end prematurely if one harmonic has no contribution. An improved criterion based on local selection of harmonics per degree-of-freedom was proposed by Grolet [START_REF] Grolet | On a new harmonic selection technique for harmonic balance method[END_REF] to overcome the limitations of previous procedures. However, this criterion requires a more complex numerical implementation. The procedure described here applies to quasi-periodic bi-dimensional HBM. For the sake of simplicity of implementation, it also relies on the approximation of the nonlinear forces but differs from [START_REF] Laxalde | Etude d'amortisseurs non-linéaires appliqués aux roues aubagées et aux systèmes multi-étages[END_REF] in the way harmonics are selected. It consists in three steps.

Once a converged solution of the nonlinear system ( 17) is obtained for a given set of harmonics with the Newton-Raphson procedure described in section 4. for instance is not retained due to a too small contribution.

The nonlinear system ( 17) is solved again using this new set of harmonics and the three steps are carried out again. This procedure is repeated until convergence of the harmonics. In general, one or two iterations are requiered for convergence.

In the exemple presented here, this automatic selection procedure converged to a set of harmonics that verifies the empirical formulae (2) and (3). This proves the reliability of the procedure. With this procedure, only harmonics that really contribute to the response are used resulting to faster and more accurate simulations. Both periodic and quasi-periodic HBM algorithms, with pseudo arc-length continuation and automatic harmonic selection, have been implemented in Code Aster via Python user scripts which give access to Aster objects and commands like FEM elementary matrix computations.

Application to the Jeffcott rotor

Data from the time transient simulations, corresponding to one point of the branch, has been used to initialise the quasiperiodic HBM. Then the pseudo arc-length continuation procedure allows following the solution branch. Fig. 17 shows the same response curve than presented in Fig. 11 with the addition of the response curve given by the quasi-periodic HBM. The branch corresponding to the quasi-periodic partial rub phenomena is perfectly reproduced. The response curve given by the quasi-periodic HBM has been ploted in Fig. 18 as a function of Ω 1 and Ω 2 . The evolution of Ω 2 against Ω 1 as yielded by the HBM has been ploted in Fig. Fig. 17 Response curve including the quasi-periodic HBM response. Close up view. HBM stable ( ), HBM unstable ( ).

Fig. 18 Response curve given by quasi-periodic HBM as a function of Ω 1 and Ω 2 . Projections of this curve onto the three plans are also ploted.

tended pseudo arc-length continuation method is very consistent with the evolution obtained from time integration. Table 2 gives the computational time required to calculate the quasi-periodic partial rub branch using traditional time integration (Runge-Kutta) and the proposed quasi-periodic HBM. The proposed method appears to be more than seven times faster than the time integration method while providing very accurate results. It appeared that the AFT procedure is the most time consuming step. Bi-dimensional FFT are also more costly compared to the FFT used for periodic HBM. The backward whip branch should also be obtained from the quasi-periodic HBM technique as well and will be considered in the next developements.

Conclusion

The capabilities of the Harmonic Balance Method (HBM), with pseudo arc-length continuation and stability assessment, have been investigated for the computation of the steadystate behaviour of rotor-to-stator contact problems. The HBM results have been compared with classical time marching solutions and analytical results [START_REF] Jiang | Determination of the global responses characteristics of a piecewise smooth dynamical system with contact[END_REF]. The application considers a simple Jeffcott rotor, with a circular, rigid, static stator. It has been demonstrated that although the rotor is submitted to unbalance only, there are four possible kinds of solution: two of these are periodic (no-rub motion and full annular rub), whereas the remaining two are quasi-periodic (partial rub and backward whip).

The periodic HBM is shown to be capable of accurately predict the two periodic solutions. The periodic HBM can also predict initiation of the quasi-periodic partial rub regime, which corresponds to a loss of stability of the periodic solution. However, the partial rub and backward whip motions cannot be described by means of the periodic HBM in the presented form.

An extension of the HBM to the treatment of quasi-periodic solutions has been proposed. This new method automati-cally determines the second fundamental frequency of the system which is the result of a self-excitation caused by contact. An extension of the pseudo arc-length continuation method made possible the continuation of a quasi-periodic branch despite the second frequency having an a priori unknown evolution along the branch. In addition, in order to improve both performance and accuracy of simulations, an automatic harmonic selection procedure has been proposed. With this procedure, only predominant harmonics are automatically used during simulation, requiring no intervention from the user. This method allowed to determine the partial rub branch and its results were very consistent with time integration results at a much lower cost.

The proposed quasi-periodic HBM currently requires to be initialised with data from time integration simulations. However, this requirement could be removed with the help of a branch switching procedure [START_REF] Beyn | Numerical continuation, and computation of normal forms[END_REF][START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF] that would allow a smooth transition between periodic full annular rub and quasi-periodic partial rub. This possibility is currently being investigated.

Although results presented here concern mainly academic examples, the method has been successfully applied to industrial rotating machinery finite-element models [START_REF] Peletan | Stratégie de modélisation simplifiée et de résolution accélérée en dynamique non linéaire des machines tournantes[END_REF].

Fig. 3

 3 Fig. 3 Response curve for µ = 0.2.

Fig. 4 Fig. 5

 45 Fig. 4 Orbit during the no rub motion at Ω = 0.14

Fig. 6 Fig. 7

 67 Fig. 6 Orbit during partial rub at Ω = 0.33

Fig. 8 Fast

 8 Fig. 8 Fast Fourier transform of the nonlinear contact force signal during partial rub at Ω = 0.33.

Fig. 9

 9 Fig. 9 Evolution of the second fundamental frequency Ω 2 as a function of Ω 1 with error bars.

Fig. 10

 10 Fig. 10 Pseudo arc-length continuation.

Fig. 11

 11 Fig. 11 Response for µ = 0.2. HBM stable ( ), HBM unstable (•).

Fig. 12

 12 Fig. 12 Real and imaginary parts of Floquet's multipliers for the unstable solution at Ω ≈ 0.33 (secondary Hopf bifurcation).

Fig. 13

 13 Fig.13Pseudo arc-length continuation with two parameters linked by an unknown function.

e 2 =

 2 [0, . . . , 0 n , 0, . . . , 0, m 12 , 0, . . . , 0 2n , . . . , 0, . . . , 0, m j2 , 0, . . . , 0 2n , . . . , 0, . . . , 0, m N2 , 0, . . . , 0 2n , 0, 0] (

Fig. 14 Fig. 15

 1415 Fig. 14 Automatic harmonic selection, step 1.

DetectionFig. 16

 16 Fig. 16 Automatic harmonic selection, step 3.

2 :- 1 :

 21 Step Bi-dimensional FFT of nonlinear forces. The nonlinear forces are calculated in the time domain from displacements. In Fig.14, the nonlinear forces are plotted as a function of θ 1 and θ 2 besides their FFT representation where each pixel represents a frequency component. A blue pixel represents a component with a very low contribution and a red pixel a high energy compo-nent. The spectrum yielded by FFT has aliasing and information is redundant.-Step 2: Suppression of redundancy (see Fig.15) and spectrum shift. The first half of rows or the first half of columns of the spectrum may be arbitrarily kept. For instance, here, the first half of rows are kept and the right half of the first row removed. For convenience, the right half of the spectrum is shifted with the left half. -Step 3: Harmonics detection (see Fig.16). The components with the highest contribution are extracted from the spectrum and added, one by one, until the error between nonlinear forces in the time domain and their reconstruction from the remaining harmonics equals a threshold value. The error is a norm of the difference between the two signals computed on the whole time vector. In other words, only important harmonics are retained. In the given example, sixteen harmonics are detected (indicated by an orange dot): {(0,1); (1,0); (1,2); (2,1); (2,3); (3,2); (3,4); . . . }. Despite being not nil, harmonic (4, 5)

19

 19

  . The evolution of Ω 2 as a function of Ω 1 found by the ex-Produit par HAL -16 Sep 2014 P a r t i a l r u b P a r t i a l r u b F u ll a n n u la r r u b F u ll a n n u la r ru b

Fig. 19

 19 Fig. 19 Evolution of the second fundamental frequency Ω 2 as a function of Ω 1 .

Table 1 Ω

 1 na

			)
	0	1	2.08
	1	4.08	5.16
	2	7.16	8.24
	3	10.2	11.3

and Ω nb calculated for the first four values of n for partial rub motion at Ω = 0.33. Here, Ω 1 = 1 × Ω and Ω 2 ≈ 2.08 × Ω.

n Ω na (×Ω) Ω nb (×Ω

Table 2

 2 Comparison of the CPU time required for the computation of the partial rub branch using the proposed method and traditional time integration.

		Time	Quasi-periodic
		integration	HBM
	Nb. of unknowns	2	67
	Nb. of points in branch	36	274
	CPU time (s)	1336	184 (7.30 speed-up)

Produit par HAL -16 Sep 2014

Acknowledgements This work was partially supported by the French National Agency (ANR) in the framework of its Technological Research COSINUS program (IRINA, project ANR 09 COSI 008 01 IRINA)


Tim Nagel
email: tim.nagel@grenoble-inp.fr

Julien Chauchat
email: julien.chauchat@grenoble-inp.fr

Achim Wirth
email: achim.wirth@legi.cnrs.fr

INVESTIGATION OF THE MULTI-SCALE INTERACTIONS BETWEEN AN OFFSHORE WIND TURBINE WAKE AND THE OCEAN-SEDIMENT DYNAMICS IN A INDEALIZED FRAMEWORK

Keywords: Offshore Wind Turbine, Wake, Ocean, Seabed, Dynamics

A coupled two dimensional idealized numerical model of the ocean and sediment layers, forced by an offshore wind turbine wake is used to investigate the complex interactions between the wake, the ocean and the sediment layers, together with the retroaction on the wind energy. Results show that the turbine wake has an impact on both, the ocean and the sediment layers. The turbine wake impacts the ocean surface and generates instabilities or vortex streets for some parameter values. Shallow ocean layers (typically below 15m) are laminar. When water depth is higher, large scale instabilities are generated, leading to a turbulent dynamic in the ocean layer. The size of the generated vortices in the ocean increases with water depth and decreases with the quadratic-law bottom friction coefficient. Considering the morphodynamics three cases are observed, depending on whether the ocean dynamics is laminar (i), has a localized (ii) or domain wide (iii) turbulent behavior. In the first case, changes in seabed elevation are around a few millimeters per month. Results are similar for the localized turbulence case with small spatial variations. For the domain wide turbulence case (iii), instantaneous seabed changes are of the order of a few millimeters per month, whereas the transport averaged over several days decreases to a few tenths of millimeter per month. This behavior is easily explained by the oscillating local velocity which transports sediments back and forth. The above emphasizes that the water depth is a key parameter for the coupled atmosphere-ocean-sediment system around wind turbines. Furthermore, considering the ocean velocity in the atmospheric forcing at the ocean surface leads to a decrease of 4 % of the power lost by friction at the atmosphere-ocean interface. Ocean dynamics could thus have a non-negligible feedback on the wind power available for the turbines and its variability.

INTRODUCTION

Because of the rising need for sustainable energy, and because up to now wind energy is one of the few forms of renewable energy that can be harvested efficiently, many European countries are planning and building offshore wind farms to increase the proportion of renewable energy in their energy mix. According to the European Wind Energy Association 2013 annual report [START_REF] Ewea | Deep water, the next step for offshore wind energy[END_REF], the installed European capacity was 5 GW at the end of 2012, and by 2020 it could be multiplied by 8 (40 GW), corresponding to 4% of the European electricity demand. By 2030, offshore wind capacity could totalize 150 GW, corresponding to 14% of the actual EU's total electricity consumption.

The environmental impact of wind farm construction is the subject of an important and diversified literature. A consequent part of this literature focus on the direct impacts of the solid structures of the wind farm on the oceanic environment and principaly on the fauna. An important environmental consequence is the birds and bats collision with the turbines [START_REF] Tucker | A mathematical model of bird collisions with wind turbine rotors[END_REF], [START_REF] Drewitt | Assessing the impacts of wind farms on birds[END_REF]). Concerning the seabed, it is affected by scour process due to the pile presence, a phenomenon similar to that of occurs at bridge piers [START_REF] Breusers | Local scour around cylindrical piers[END_REF], [START_REF] Breusers | Scouring[END_REF], [START_REF] Roulund | Numerical and experimental investigation of flow and scour around a circular pile[END_REF]). It has also been shown that a local seabed elevation occurs in offshore wind farms [START_REF] Van Der Veen | Seabed morphodynamics due to offshore wind farms[END_REF]. On a larger scale, impact of wind farms on the European regional climate [START_REF] Vautard | Regional climate model simulations indicate limited climatic impacts by operational and planned european wind farms[END_REF] or on hurricanes [START_REF] Jacobson | Taming hurricanes with arrays of offshore wind turbines[END_REF] have been studied recently, showing that the wind farms environmental impacts are an important questions nowadays.

Nevertheless, most of the interactions between offshore wind farms and their local environment need further research and the understanding of the global impact induced by the construction of offshore wind farms in coastals areas may requires the modeling of the ocean-atmosphere-sediment coupled system. However, if air-sea and ocean-sediment interactions have been the subject of an important literature, to the best of our knowledge, no study has been done on the multi-scale interactions in an atmosphereocean-sediment coupled system. This is even more true in the vicinity of Marine Renewable Energy Devices. This paper shows the first results of an idealized 2D numerical model build to study the impact of an offshore turbine wake on the ocean and the sediment dynamics. In Section 2., the physical model and the basic equations that describe the mathematical model are detailled. The numerical model used here is described in Section 3., while Section 4. is devoted to the presentation and the discussion of the results. The summary and conclusions are given in Section 5.

PHYSICAL AND MATHEMATICAL MODEL

The physical model consists in two superposed layers (figure 1.a), an homogeneous shallow water ocean layer above a sediment bed layer, composed of cohesion-less particles. The atmospheric layer is represented as an external forcing (F ), which corresponds to the wake of wind turbines. The domain length in the x-direction is denoted as Lx and as Ly in the y-direction. The average depth of the ocean layer is denoted as H. The local thickness of the ocean layer and the seabed elevation, are denoted as h(x, y, t) and hs(x, y, t) respectively. The free surface elevation is:

η(x, y, t) = h(x, y, t) -H + hs(x, y, t) [1]
The dimensional variables are consigned in the table 1. Densities are denoted as ρ, ρs and ρa for the ocean, the sediment and the atmosphere, respectively, Ds is the diameter of a sand grain. The ocean layer is forced by the local wind stress at its upper surface.

The spatial variation of this wind stress incorporates the wake-profile of a wind turbine. The oceanic motion induces a shear stress τ b on the sediment bed layer. This stress and the seabed elevation are responsible for the coupling between the ocean and the sediment bed layers.

Hydrodynamic model

Mathematically, the ocean dynamic is given by the two dimensional Shallow Water (SW) equations, as most of the wind farms are localized in coastal areas (in 2012 the average water depth of offshore wind farms was 22m and the average distance to shore was 29km (EWEA, 2013)), such equations seem to be particularly adapted to the problem:

∂t -→ u + ( -→ u • ∇) -→ u + g∇η = ν∇ 2-→ u + -→ F - -→ G [2] ∂th + ∇(h -→ u ) = 0, [3]
where -→ F and -→ G correspond to the frictionnal accelerations applied on the ocean layer by the atmosphere and the sediment, respectively. These frictionnal accelerations correspond to the forces excerted on the ocean by the wind shear and the bottom shear respectively and are characterized by a quadratic friction law. The atmosphere-ocean friction law ( -→ F ) is defined by:

-→ F = 1 ρh -→ f , [4]
where -→ f is the shear stress applied to the ocean. The latter is calculated using the velocity difference between wind and ocean current as described by [START_REF] Moulin | A drag-induced barotropic instability in air-sea interaction[END_REF]:

-→ f = CD ρa || -→ u a --→ u o || ( -→ u a --→ u o) [5] with || -→ u a --→ u o ||= (ua -uo) 2 + (va -vo) 2 .
The drag coefficient CD follows [START_REF] Wu | Wind-stress coefficients over sea surface from breeze to hurricane[END_REF] and [START_REF] Smith | Coefficients for sea surface wind stress, heat flux, and wind profiles as a function of wind speed and temperature[END_REF] suggestion: CD = (0.6 + 0.07u).10 -3 for 6 m/s < ua < 26 m/s

[6]

The validity range of this friction law corresponds to the good working conditions of wind turbines (from 4 m.s -1 to 25 m.s -1 ). A similar quadratic friction law is used to model the friction between the ocean and the sediment layers:

-→ G = 1 h -→ τ b , [7] -→ τ b = CDs || -→ u || -→ u , [8]
with || -→ u ||= √ u 2 + v 2 . CDs = 0.005 correspond to the friction coefficient between the sand and the ocean.

Morphodynamic model

For a bed composed of cohesionless grains (equivalent to sand), the sediment starts to move when the drag force exerted by the flow is higher than the friction force between the grains. The dimensionless Shields parameter [START_REF] Shields | Anwendung der ahnlichkeits mechanik und der turbulenzforschung auf die geschiebebewegung[END_REF], gives the ratio between the drag force and the apparent submerged grains weight:

θ = τ b (ρs -ρ)gd , [9]
where g is the gravitational acceleration and d is the characteristic sediment particle diameter.

The sediment starts to move as soon as the Shields parameter exceeds a typical critical value (θ > θc). The critical Shields parameter, θc gives the threshold of motion for sediments at the bed and depends on several sediment properties, as the density or the grain size.

Bed motion results from a local flux balance. The mass conservation equation, also called Exner equation, allows to calculate the time evolution of the seabed elevation (hs):

∂ths(x, y, t) + -→ ∇ • -→ q (x, y, t) = 0, [10]
where -→ q is the total sediment flux.

Two modes of sediment transport exists, bed-load transport and suspended load transport. In this work, only the bed-load transport will be considered, this type of transport is generally dominating for rather low bed shear stress, i.e when the Shields number of the flow is juste above the critical value.

If only bed-load transport is considered, the total sediment flux becomes: [11] where p is the bed porosity and -→ q b is the bed load transport rate.

-→ q = 1 1 -p -→ q b ,
In the present model, the Meyer-Peter and Müller (1948) (MPM) transport formula is used. It relates the particle flux q b to the excess Shields parameter:

q b d (ρs/ρ)gd = 8(θ -θc) 3/2 if θ > θc 0 otherwise [12]
2.3 Turbine wake model According to [START_REF] Barbe | Etude de l'impact des éoliennes offshores sur l'ocean et l[END_REF], the turbine wake affects the ocean by decreasing the surface shear in the region where the wake intersects the ocean surface (see figure 1). This perturbation is maximal at the impact location and then decreases downstream. A distance of 40 to 60 times the rotor diameter D is necessary for the wake effects to vanish.

The wake model used in this work describes the velocity deficit induced by the rotor in the far wake region. The simplest model is the [START_REF] Jensen | A note on wind generator interaction[END_REF] one, assuming a linearly expanding or cone wake with a velocity deficit which only depends of the distance to the rotor.

u = U∞ 1 - 1 - √ 1 -Cw (1 + 2kx D ) 2 , [ 13
]
where U∞ is the wind velocity far from the turbine, Cw the drag coefficient between the turbine and the air and k the Wake Decay Constant. For offshore turbines, the admitted value is equal to 0.05.

The wake impacts the ocean surface at a given distance downstream from the rotor position (figure 1.b). For the Jensen's model we have:

ximpact = Hr -D/2 k [14]
This impact distance depends on the wind turbine height and on the rotor diameter and is increasing with the turbine size. For the turbines considered here (Hr = 70m and D = 80m, k = 0.05), the impact distance is equal to 600 meters, several times the turbine height. Perturbations in the seabed induced by the turbine pile are localized in the pile vicinity, the seabed is thus perturbed on a local scale. As the pile diameter doesn't exceed 5 meters for the type of turbine modeled, seabed perturbations induced by the pile and the wake presence are uncorrelated. The effect of the pile is thus not considered here.

The Jensen wake model's boundaries are extremely sharp, indeed, the velocity difference between inside and outside the wake corresponds to a step function. Such sharp boundaries are unrealistic and are prone to generating instabilities. To be more realistic without changing the large scale characteristic of the wake model, a good solution consists in using a mollifier (also known as approximation to the identity) function. Thus, close to the wake boundary, the wind velocity ua is smoothed and becomes ua by using the following operator:

ua(y) = 1 Lm √ π b a ua(y )e -y-y Lm 2 dy , [ 15
]
where a = Ly/2 and b = Ly or b = 1 are the integral bounds, for the upper and the lower wake boundaries, respectively, y = y (x) is the position of the wake boundary and Lm ≈ 1 -10 dy is the mollify length, for Lm = 1, the transition zone correspond to a step function.

NUMERICAL MODEL

The ocean and the sediment are considered in a rectangle of size Lx × Ly. The numerical grid is regular and contains nx × ny points. The number of grid points is equal to the length (in meters) of the domain in each direction, i.e Lx = 2000m and Ly = 600m.

The spatial resolution is the same in both horizontal directions, ∆x = Lx/nx = ∆y = Ly/ny = 1 m. Periodic boundary conditions are used in both horizontal directions.

The overall model can be divided in two coupled modules, the hydrodynamic and the morphodynamic ones. Concerning the hydrodynamic module, where the shallow water equations for the ocean are solved, a second-order, centered, finite-difference method is used for the space discretization, and a second-order Runge-Kutta scheme is used for the time discretization. The morphodynamic module solves the Exner equation using a NOCS (Non-Oscillatory Central Scheme) scheme as described by [START_REF] Jiang | High-resolution nonoscillatory central schemes with nonstaggered grids for hyperbolic conservation laws[END_REF].

For a given spatial resolution, the time step is imposed by the Courant-Friedrichs-Lewy (CFL) condition. The sand waves propagation celerity corresponds to dunes migration speed which is approximatively 1000 times slower than oceanic waves. It is thus the CFL condition for the ocean that sets the maximum time step, here ∆t = 0.02 s. As the morphological timescale is large compared to the hydrodynamic one, two time-steps are used, the morphological the time-step is ∆t morpho = 1000 ∆t for all the simulations.

RESULTS AND DISCUSSIONS

All the numerical simulations presented in this work involve the same wake parameters, such as the rotor diameter or the wind velocity. The rotor diameter and height (D = 80m, Hr = 70m) correspond to a Vestas V80-2.0 MW wind turbine, one of the most widespread offshore wind turbines. The wind speed is equal to 20m.s -1 , corresponding to the high range part of turbine good working conditions (from 4m.s -1 to 25m.s -1 ) 1 . Simulations have been undertaken for five different water depths and a summary of the different simulations parameters is given in table 2. Run H20CD2P3 has been undertaken with a bottom friction coefficient C Ds = (2/3) * CDs, run H15WUo without the ocean velocity in the calculations of the atmospheric forcing, run H30WUow without the wake and the ocean velocity in the calculations of the atmospheric forcing and run H30Ww without including the wake in the calculations of the atmospheric forcing. So in run H30WUow and H30Ww the forcing is constant in time and space and the solutions in the ocean and sediment converge to a stationary state (and can be calculated analytically). For each case, the initial free surface and seabed elevation fields are equal to zero.

(a) (b) (c) (d) (e) (f) (g) (h) (i) (j)

Generation of instabilities and impact on morphological evolution

The results presented in figure 2 show that increasing water depth leads to the generation of instabilities in the wake. The ocean dynamics becomes turbulent, starting from a vortex formation at the wake boundaries. Furthermore, when eddies are generated, their size and spacing depend also on the water depth, as summarize in table 2.

For the different water depth studied here, three cases of ocean dynamics are observed. Up to 15 meters of water depth, the dynamics is laminar, the vorticity is higher at the wake boundaries than in other parts of the domain but no instabilities are generated (figure 2.a). When the water depth is between 20 and 30 meters, vortices are generated at all the wake boundaries (see figures 2.c and 2.e). Two independents localized vortex streets appear. The third case emerges when water depth is higher than 40 meters.

As the size of the vortices increases, the two vortex streets interact leading to a domain wide turbulence (see figures 2.g and 2.i).

An increase of approximately a factor 2 of water depth leads to a change of the same order in the vortex spacing (Table 2).

Decreasing the bottom friction coefficient in order to conserve the friction coefficient and water depth ratio (S = CD/h) for two different water depth (cases H20CD2P3 and H30) gives similar results (figure 3.a and figure 3.b). It appears thus that vortex size and spacing seems to depend on S . This ratio, S, appears in the bottom friction terms in the SW equations:

-→ G = 1 h -→ τ b = S || -→ u || -→ u [16]
This result confirms the findings of Chen and Jirka (1997) on the cylinder wake in a SW flow where S = (D * CD)/h, with D the diameter of the cylinder. In the present paper no dependency test upon the wake width at the impact location (corresponding to the diameter of the cylinder of Chen and Jirka (1997) work) have been undertaken. Here, for increasing water depth, S is decreasing and the bottom friction force as well, allowing for stronger instabilities. For all the simulations undertaken, the qualitative spatial long-time impact on the wake upon the seabed is similar (figure 2, right column). Figure 2 shows that the velocity deficit induced by the wake leads to a sediment accumulation, i.e a dune formation close to the impact (between x = 300m and x = 500m) which corresponds to the local wake width. The same phenomenon occurs at the wake boundaries. Such formations are due to the non uniformity of the oceanic velocity field and thus of the local bottom shear stress. Downstream in the wake, the velocity deficit becomes less important, corresponding to a bottom shear stress that induce bed erosion. Outside the wake, the flow velocity is higher, increasing thus the bottom shear stress and thus the sediment transport, leading to a hole formation on each side of the dune. After x = 1000m, the velocity and the bottom shear stress decrease along the x-direction, stopping the erosion and starting sediment accumulation when moving downstream in this region (figure 2, right column).

Figure 4 shows the bathymetry variation between two consecutive outputs (∆t = 2000 s) of the same simulation. The spatial pattern of the morphodynamical evolution between the turbulent and the laminar cases is completely different. The vortices strongly affect the short term seabed morphodynamics, as expected. For cases H15 and H20 (figures 4.a and 4.b, respectively), the bathymetry variations between two consecutive outputs are qualitatively similar to the seabed elevation fields after 11.5 days of dynamics. This is not the case from 30 m to 50 m water depth (figures 4.c and 4.d, respectively), the wakes imprint is totally annihilated by the signature of the large scale vortices. Difference in vortex scale observed in the ocean is recovered in the seabed, showing that vortex formation in the ocean has an significant impact on the seabed morphodynamics.

Considering the morphodynamical evolution, three cases can be observed, depending on whether the ocean dynamics is laminar, has a localized or a domain wide turbulent behavior. In the first case, changes in seabed elevation are around a few millimeters per month (figure 2.b). Results are similar for the localized turbulence case with small spatial variations (figure 2.d). For the domain wide turbulence case, instantaneous seabed changes are of the order of a few millimeters per month (figure 4.d), whereas the transport averaged over several days decreases to a few tenths of millimeter per month (figure 2.j). This behavior is easily explained by the oscillating local velocity which transports sediments back and forth. This result may cause environmental issues. Sediment motion induced by the turbulent state in the ocean could increase the local turbidity and perturb the photosynthesis in the water column even when the long term accumulated transport is small.

Air-sea interactions and energetic balance

In this idealized model the source of mechanical energy is the atmospheric dynamics. Different computations have been undertaken using different water depth in order to quantify the energy transfer, the results are presented in figure 5. In this analysis, we neglected the y component of the velocity vectors, they are vanishing for the atmosphere and subdominant for the ocean.

The power lost by the atmosphere (Pa) is the integral over the domain of the product between the wind stress and the wind velocity.

Pa = CDρa S || -→ u a --→ u o || ( -→ ua --→ uo) • -→ ua dS.
[17]

As the y component of the velocity is vanishing in the atmosphere and subdominant in the ocean (vo uo, ua), equation ( 17) simplifies to:

Pa = CDρa S | ua -uo | (ua -uo)ua dS.
[18] Furthermore, as the velocity in the atmosphere is always larger than in the ocean (ua > uo, ∀x, ∀y, ∀t):

Pa = CDρa S (ua -uo) 2 ua dS.
[19]

The power taken by the ocean (Po) is the integral over the domain of the product between the wind stress and the ocean velocity, using the above simplifications it can be written as:

Po = CDρa S (ua -uo) 2 uo dS.
[20]

According to the Betz theory, the maximal power value which can be taken from the wind by the turbine is lower or equal to 16/27 of the wind power without the turbine presence: [21] where R correspond to the radius of the rotor. Figure 5 shows that for all the cases investigated, the power received by the ocean from the atmosphere is several orders of magnitude lower than the power lost by the atmosphere (figures 5.a and 5.c). The major part of the power lost by the atmosphere is thus dissipated by the friction between the atmosphere and the ocean. Taking into account the ocean velocity in the wind forcing leads to a decrease by 2% of the power received by the ocean from the atmosphere (comparison between cases H30WUow, H30Ww and cases H15, H15WUo of figure 5.c) and by 4% of the power lost by the atmosphere (comparison between cases H30WUow H30Ww and cases H15, H15WUo of figure 5.a). For case H15, the loss of atmospheric power due to the friction is reduced by approximatively 4% as compared to case H15WUo. Taking into account the relative velocity between the atmosphere and the ocean leads to a reduction of the friction between these two layers and thus to the power loss by the atmosphere. Furthermore, considering the atmospheric wake presence leads to a decrease of the power received by the ocean and the one lost by the atmosphere of approximatively 16%. Such phenomenon may be explain by the fact that the wake impact reduces the local atmospheric velocity more than the local ocean one, the relative velocity between the two layers is thus decreasing, so as the friction. For all the depths considered, the power loss to friction, which is the difference between the power lost by the atmosphere and the one received by the ocean, is higher than the power taken by the wind turbine. The difference between the values for the laminar (H15) and the domain wide turbulence case (H50) is small, of the order of 0.01% for the power lost by the atmosphere and of the order of 0.3% for the power received by the ocean. It shows that the oceanic turbulence as subdominant role in the air-sea energetic balance.

P turbine = 16 27 1 2 ρaπR 2 (ua) 3 ,
An analytical development can confirm and explain the above observations. Time average and developed equation 19 can be written for laminar and turbulent oceanic flow, respectively:

Pa l = CDρa S ua 3 -2ua 2 uo + uauo 2 dS.
[22]

Pa t = CDρa S ua 3 a 1 -2ua 2 uo + uauo 2 a 2 + uau o 2 a 3 dS, [ 23
]
where uo is the mean flow velocity and u o 2 is the square of the velocity fluctuations, or the x-normal component of the Reynolds stress tensor.

The same procedure can be applied to the power received by the ocean (20) and gives for laminar and turbulent oceanic flow, respectively:

Po l = CDρa S ua 2 uo -2uauo 2 + uo 3 dS.

[24]

Po t = CDρa S ua 2 uo o 1 -2uauo 2 + uo 3 o 2 + u o 2 (3uo -2ua) o 3
dS.

[25]

The importance of the different terms in the air-sea energetic balance appears in the equations ( 22), ( 23), ( 24) and ( 25). The terms a1 = ua 3 and o1 = ua 2 uo are linked to the power taken by the ocean without taking into account the oceanic velocity in the atmospheric forcing. These terms are predominant in laminar and turbulent equations. The terms a2 = -2ua 2 uo 2 + uauo 2 and o2 = -2uauo 2 + uo 3 correspond to mean flow power correction occurring when the oceanic velocity is taking into account in the atmospheric forcing. As the velocity in the atmosphere is always much larger than in the ocean, these correction terms are negative. Taking into account the oceanic velocity in the atmospheric forcing leads thus to reduce both Pa and Po.

For equations ( 23) and ( 25) the terms a3 = uau o 2 and o3 = u o 2 (3uo -2ua) correspond to the contributions of the turbulent fluctuations in the energetic balance. The terms a3 > 0 and o3 < 0, at least for the idealized situation considered here. When turbulence occurs in the domain, it increases thus the power lost by the atmosphere but decreases the power received by the ocean, a result in agreement with the observations done in figures 5.b and 5.d. This can be explain by the fact that a part of the turbulent fluctuations is not in the same direction as the oceanic motion and the atmospheric forcing and increases thus the friction between the two layers. Furthermore, these fluctuations certainly reduce the energetic transfer efficiency to the mean flow which may explain that they decrease the power received by the ocean.

If the turbulence seems to play a minor role in the air-sea energetic balance, the results presented in Section 4.1 show its importance on both oceanic and morphological dynamics. The equations ( 23) and ( 25 Figure 6 shows the 2D Urms fields for different cases investigated. For each case, the maximum of Urms is localized at the wake boundaries, where eddies are generated. This is particularly clear for case H20 and H20CD2P3 (figures 6.a and 6.b, respectively)

where the vortex streets are independent. When the turbulence is domain wide (figure 6.c), the highvalues of Urms are less localized and start to fill large parts of the domaine. Furthermore, the maximum value of the Urms is increasing for increasing water depth. The maximum value of the Urms is equal to 0.00957 m.s -1 (for the H50 case), which is only 2.4% of the mean velocity (of the order of 0,4 m.s -1 ), showing once again the fluctuation's small order of magnitude compared to the mean flow.

Computations of the Vrms = v o 2 velocity (figure 7) show that the y-component of the turbulent fluctuations is of the same order of magnitude than the x ones. Variations of free surface elevation for the different computations undertaken are extremely small, of the order of thenths of a millimeter. The dynamic height is thus small, a result in agreement with potential energy calculations showing that potential energy involved in the problem is eight orders of magnitude lower than kinetic energy. The major part of the energy involved in the problem is thus kinetic and not potential.

SUMMARY AND CONCLUSION

An idealized 2D numerical model has been proposed to study the impact of an offshore wind turbine wake on the ocean and sediment dynamics. To the best of our knowledges so far no study has been done on this subject. A simplified physical model has been proposed and a mathematical model has been consistently built.

The present results show that the turbine wake has an impact on both ocean and sediment bed layers. Turbine wake impact on the ocean surface can generate instabilities and vortex streets formation. Size and spacing between these vortices seems to be controlled by the friction between the bed and the ocean layer and more precisely by the friction coefficient and water depth ratio, S = CD/h. When S is decreased, large scale instabilities are more easily generated, leading to a domain wide turbulence state in the ocean. Concerning the seabed, results show that the non-uniformity of the oceanic velocity field induced by the wake presence is recovered in the local bottom shear stress responsible for sand erosion and deposition. As the vortices strongly affect the seabed morphodynamics, for localized and domain wide turbulence, the wakes imprint tends to be reduced by large scale vortices and the oscillating local velocity which transports sediments back and forth.

Taking into account the ocean velocity in the wind forcing leads to a decrease by approximatively 4% of the power loss to friction of the atmosphere, for the case considered here, by reducing the relative velocity between the two layers. Power received by the ocean is also decreasing by approximatively 2%. Finally, considering the wake presence leads to a decrease of 16% of the power received by the ocean, as the local air-sea friction is decreased.

The results also show that even if the turbulence is strongly influencing both, oceanic and morphological dynamics its role is minor in the air-sea energetic balance. Furthermore, as the atmospheric dynamics is not resolved here, the oceanic turbulence has a limited retro-action on the forcing, considering the atmosphere as a shallow-water layer may increase the importance of these interactions.

Figure 1 :

 1 Figure 1: Sketch of the physical model, exponent "s" refers to the sediment bed (a) and of the Jensen wake (b).

Figure 2 :

 2 Figure 2: 2D vorticity (s -1 , left column) and seabed elevation (right column) fields after 11.5 days of dynamics. Cases H15 (a and b), H20 (c and d), H30 (e and f), H40 (g and h) and H50 (i and j).

Figure 3 :

 3 Figure 3: 2D vorticity seabed elevation fields after 11.5 days of dynamics. Cases H30 (a) and H20CD2P3 (b).

Figure 4 :

 4 Figure 4: 2D bathymetry difference field between two consecutive outputs for H15 (a), H20 (b), H30 (c) and H50 (d) after 11.5 days.

Figure 5 :

 5 Figure 5: Power lost by the atmosphere for the different cases investigated and by the wind turbine (a), for different turbulent states (b) and power received by the ocean for the different cases investigated (c) and for different turbulent states (d).

  Figure5shows that for all the cases investigated, the power received by the ocean from the atmosphere is several orders of magnitude lower than the power lost by the atmosphere (figures 5.a and 5.c). The major part of the power lost by the atmosphere is thus dissipated by the friction between the atmosphere and the ocean. Taking into account the ocean velocity in the wind forcing leads to a decrease by 2% of the power received by the ocean from the atmosphere (comparison between cases H30WUow, H30Ww and cases H15, H15WUo of figure 5.c) and by 4% of the power lost by the atmosphere (comparison between cases H30WUow H30Ww and cases H15, H15WUo of figure 5.a). For case H15, the loss of atmospheric power due to the friction is reduced by approximatively 4% as compared to case H15WUo. Taking into account the relative velocity between the atmosphere and the ocean leads to a reduction of the friction between these two layers and thus to the power loss by the atmosphere. Furthermore, considering the atmospheric wake presence leads to a decrease of the power received by the ocean and the one lost by the atmosphere of approximatively 16%. Such phenomenon may be explain by the fact that the wake impact reduces the local atmospheric velocity more than the local ocean one, the relative velocity between the two layers is thus decreasing, so as the friction. For all the depths considered, the power loss to friction, which is the difference between the power lost by the atmosphere and the one received by the ocean, is higher than the power taken by the wind turbine. The figures 5.b and 5.d show the power lost by the atmosphere and received by the ocean, respectively, for runs H15, H20, H30, H40 and H50. Values of the power lost are increasing for increasing water depth and thus increasing turbulence level (figure 5.b) but the values of power received by the ocean are decreasing for increasing water depths (figure 5.d). The difference between the values for the laminar (H15) and the domain wide turbulence case (H50) is small, of the order of 0.01% for the power lost by the atmosphere and of the order of 0.3% for the power received by the ocean. It shows that the oceanic turbulence as subdominant role in the air-sea energetic balance.

  ) seem to show that evaluating the x-component of the root-mean-square velocity (Urms = u o 2 ) of the fluctuations can be used to estimate the importance of turbulent fluctuations in the ocean for the energy fluxes at the atmosphere-ocean interface.

Figure 6 :

 6 Figure 6: 2D Urms fields for H20 (a), H20CD2P3 (b) and H50 (c) after 11.5 days of dynamics.

Figure 7 :

 7 Figure 7: 2D Vrms fields for H20 (a), H20CD2P3 (b) and H50 (c) after 11.5 days of dynamics.

Table 1 :

 1 Domain parameters.

	Ly (m)	H (m)	Ds (m)

ρ (kg.m -3 ) ρs (kg.m -3 ) ρa (kg.m -3 ) Lx (m)

Table 2 :

 2 Wake simulations and their parameters.

	Run name	H (m)	Eddies Diameter (De (m))	Eddies spacing (Le (m))
	H15	15	-	-
	H20	20	100	400
	H30	30	250	500
	H40	40	300-400	670
	H50	50	550	1000
	H20CD2P3	20	250	500
	H15WUo	15	-	-
	H30WUow	30	-	-
	H30Ww	30	-	-

Vestas V80-2.0 MW product brochure: http://www.vestas.com

Acknowledgements

We thank Cyrille Bonamy for the parallelization of the code and for his precious help in the numerical implementation and the Institut Carnot Énergies du Futur for their financial support. All the computations presented in this paper were performed using the Froggy platform of the CIMENT infrastructure (https://ciment.ujf-grenoble.fr), which is supported by the Rhône-Alpes region (GRANT CPER07-13 CIRA) and the Equip@Meso project (reference ANR-10-EQPX-29-01) of the programme

Investissements d'Avenir supervised by the Agence Nationale pour la Recherche.


Marco Corneli

Pierre Latouche

Fabrice Rossi

Fabrice Rossi Modelling

Modelling time evolving interactions in networks through a non stationary extension of stochastic block models

scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Modelling time evolving interactions in networks through a non stationary extension of stochastic block models

Marco Corneli, Pierre Latouche and Fabrice Rossi

Université Paris 1 Panthéon-Sorbonne -Laboratoire SAMM 90 rue de Tolbiac, F-75634 Paris Cedex 13 -France

Résumé. Le modèle à blocs stochastiques (SBM) décrit les interactions entre les sommets d'un graphe selon une approche probabiliste, basée sur des classes latentes. SBM fait l'hypothèse implicite que le graphe est stationnaire. Par conséquent, les interactions entre deux classes sont supposées avoir la même intensité pendant toute la période d'activité. Pour relaxer l'hypothèse de stationnarité, nous proposons une partition de l'horizon temporel en sous intervalles disjoints, chacun de même longueur. Ensuite, nous proposons une extension de SBM qui nous permet de classer en même temps les sommets du graphe et les intervalles de temps où les interactions ont lieu. Le nombre de classes latentes (K pour les sommets, D pour les intervalles de temps) est enfin obtenu à travers la maximisation de la vraisemblance intégrée des données complétées (ICL exacte). Après avoir testé le modèle sur des données simulées, nous traitons un cas réel. Pendant une journée, les interactions parmi les participants de la conférence HCM Hypertext (Turin, 29 Juin -1er Juillet 2009) ont été traitées. Notre méthodologie nous a permis d'obtenir une classifications intéressante des 24 heures: les moments de rencontre tels que les pauses café ou buffets ont bien été détectés. La complexité de l'algorithme de recherche, linéaire en fonction du nombre initial de clusters (K max et D max respectivement), nous oriente vers l'utilisation d'instruments avancés de classification, pour réduire le nombre attendu de classes latentes et ainsi pouvoir utiliser le modèle pour des réseaux de grand dimension.

Mots-clés. Graphes aléatoires, classification temporelle, modèles à blocs stochastiques, vraisemblance classifiante intégrée Abstract. In this paper, we focus on the stochastic block model (SBM), a probabilistic tool describing interactions between nodes of a network using latent clusters. The SBM assumes that the network has a stationary structure, in which connections of time varying intensity are not taken into account. In other words, interactions between two groups are forced to have the same features during the whole observation time. To overcome this limitation, we propose a partition of the whole time horizon, in which interactions are observed, and develop a non stationary extension of the SBM, allowing to simultaneously cluster the nodes in a network along with fixed time intervals in which the interactions take place. The number of clusters (K for nodes, D for time intervals) as well as the class memberships are finally obtained through maximizing the complete-data integrated likelihood by means of a greedy search approach. After showing that the model works properly with simulated data, we focus on a real data set. We thus consider the three days ACM Hypertext conference held in Turin, June 29th -July 1st 2009. Proximity interactions between attendees during the first day are modelled and an interesting clustering of the daily hours is finally obtained, with times of social gathering (e.g. coffee breaks) recovered by the approach. Applications to large networks are limited due to the computational complexity of the greedy search which is dominated by the number K max and D max of clusters used in the initialization. Therefore, advanced clustering tools are considered to reduce the number of clusters expected in the data, making the greedy search applicable to large networks.

Keywords. Random graphs, time event clustering, stochastic block models, integrated classification likelihood.

Introduction

Since the interactions between nodes of a network generally have a time varying intensity, the network has a non trivial time structure that we wish to infer. An example of this complexity can be observed in Figure [START_REF] Randriamanamihaga | Clustering the Vélib' dynamic Origin/Destination flows using a family of Poisson mixture models[END_REF]. On the vertical axis the aggregated number of proximity face-to-face interactions (less than 1.5 meter) between attendees of the HCM Hypertext conference (Turin, June 29th -July 1st, 2009) is given. On the horizontal axis, a time line is reported, corresponding to the 24 hours of the first day of conference. The day was partitioned in small time intervals of 20 seconds in the original data frame. We considered 15 minutes time aggregations leading to a partition of the day made of 96 consecutive quarter-hours: each red point in the figure corresponds to one of them. If we associate each attendee to a node and each interaction to an edge, the non stationarity of such a graph is clear. In the following section, after describing the general stochastic block model, we illustrate the temporal evolution we propose.

A non stationary stochastic block model

We present here the SBM (Holland et al. [START_REF] Holland | Stochastic blockmodels: first steps[END_REF]), using the same notations as in Wyse, Friel, Latouche (2014). The set of nodes A = {a 1 , . . . , a N } is introduced. Undirected links between nodes i and j from A are counted by the observed variable X ij , being the component (i, j) of the N × N adjacency matrix X = {X ij } i≤N,j≤N . Nodes in A are clustered in K disjointed subgroups respectively:

A = ∪ k≤K A k , A l ∩ A g = ∅, ∀l = g
Nodes in the same cluster in A have linking attributes of the same nature. We introduce an hidden vector c = {c 1 , . . . , c N } labelling each node's membership:

c i = k iff a i ∈ A k , ∀k ≤ K.
In order to introduce a temporal dimension, consider now a sequence of equally spaced, adjacent time steps {∆ u := t u -t u-1 } u≤U over the interval [0, T ] and a partition C 1 , . . . , C D of the same interval1 . We introduce furthermore a random vector y = {y u } u≤U , such that y u = d if and only if I u :=]t u-1 , t u ] ∈ C d , ∀d ≤ D. We attach to y a multinomial distribution:

p(y|β, D) = d≤D β |C d | d ,
where

|C d | = #{I u : I u ∈ C d }.
Now we define N Iu ij as the number of observed connections between a i and a j , in the time interval I u and we make the following crucial assumption:

p(N Iu ij |c i = k, c j = g, y u = d) follows a Poisson(∆ u λ kgd ), (1)
hence the number of interactions is conditionally distributed like a Poisson random variable with parameter depending on k, g, d (∆ u is constant).

Notation: In the following, for seek of simplicity, we will note: The adjacency matrix, noted N ∆ , has three dimensions (N × N × U ) and its observed likelihood can be computed explicitly:

p(N ∆ |Λ, c, y, K, D) = k,g,d ∆ S kgd c i c j yu N Iu ij ! e -∆λ kgd R kgd λ S kgd kgd , (2)
where we noted

S kgd := c i c j yu N Iu ij and R kgd := |A k ||A g ||C d | 2 .
The subscript u was removed from ∆ u to emphasize that time steps are equally spaced for every u.

Since c and y are not known, a multinomial factorizing probability density p(c, y|Φ, K, D), depending on a hyperparameter Φ, is introduced. The joint distribution of labels looks finally as follows:

p(c, y|Φ, K, D) = k≤K ω |A k | k d≤D β |C d | d , (3)
where Φ = {ω, β}.

Exact ICL for non stationary SBM

The integrated classification criterion (ICL) was introduced as a model selection criterion in the context of Gaussian mixture models by Biernacky et al. [START_REF] Biernacky | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF]. Côme and Latouche [START_REF] Côme | Model selection and clustering in stochastic block models with the exact integrated complete data likelihood[END_REF] proposed an exact version of the ICL based on a Bayesian approach for the stochastic block model and this is the approach we follow here. The quantity we focus on is the complete data log-likelihood, integrated with respect to the model parameters Φ and Λ = {λ kgd } k≤K,g≤K,d≤D :

ICL = log p(N ∆ , c, y, Λ, Φ|K, D)dΛdΦ . (4)
Introducing a prior distribution ν(Λ, Φ|K, D) over the pair Φ, Λ and thanks to ad hoc independence assumptions, the ICL can be rewritten as follows:

ICL = log ν(N ∆ |c, y, K, D) + log (ν(c, y|K, D)) . ( 5
)
The choice of prior distributions over the model parameters is crucial to have an explicit form of the ICL.

A priori distributions

We consider the conjugate prior distributions. Thus we impose a Gamma a priori over Λ: It can be proven that the two terms in ( 5), reduce to:

ν(λ kgd |a, b) = b a Γ(a) λ a-1 kgd e -bλ kgd
ν(N ∆ |c, y, K, D) = k,g,d b a ∆ S kgd Γ(a) c i c j yu N Iu ij ! (6)
Γ(S kgd + a) [∆R kgd + b] S kgd +a and:

ν(c, y|K, D) = Γ(αK) Γ(α) K k≤K Γ(|A k | + α) Γ(N + αK) × Γ(γD) Γ(γ) D d≤D Γ(|C d | + γ) Γ(U + γD) . ( 7
)

ICL maximization and experiments

In order to maximize the integrated complete likelihood (ICL) in equation ( 5) with respect to the four unknowns c, y, K, D, we rely on a greedy search over labels and the number of nodes and time clusters. This approach is described in Wyse, Frial and Latouche [START_REF] Wyse | Inferring structure in bipartite networks using the latent block model and exact ICL[END_REF] for a stationary latent block model.

Experiments

Simulated data

Some experiments on simulated data have initially been conducted. Based on the model described above, we simulated interactions between 50 nodes, clustered in three groups (K = 3). Interactions take place into 24 time intervals of unitary length (ideally one hour), clustered into three groups too (D = 3). Nodes and time intervals labels are sampled from multinomial distributions. With these settings, we consider 27 different Poisson parameters (λs) generating connections between nodes. The generative model used to produce them is described by:

λ kgl = s 1 [k] + s 2 [g] + s 3 [l], k, g, l ∈ {1, 2, 3}
where:

s 1 = [0, 2, 4] s 2 = [0.5, 1, 1.5] s 3 = [0.5, 1 , 1.5]
ans s 1 [k] denotes the k-th component of s 1 . Similarly for s 2 and s 3 . The greedy search algorithm we coded was able to exactly recover these initial settings, converging to the true ICL of -122410. Other experiments were run with different values inside vectors s 1 , s 2 , s 3 . Not surprisingly the more nuanced differences between λs are, the more difficult it is for the algorithm to converge to the true value of the ICL3 .

Real Data

The dataset we used was collected during the ACM Hypertext conference held in Turin, June 29th -July 1st, 2009. Conference attendees volunteered to to wear radio badges that monitored their face-to-face proximity. The dataset represents the dynamical network of face-to-face proximity of 113 conference attendees over about 2.5 days 4 . Further details can be found in Isella, Stehlé, Barrat, Cattuto, Pinton, Van den Broeck [START_REF] Isella | What's in a crowd? Analysis of face-to-face behavioral networks[END_REF]. We focused on the first conference day, namely the twenty four hours going from 8am of June 29th to 7.59am of June 30th. The day was partitioned in small time intervals of 20 seconds in the original data frame. We considered 15 minutes time aggregations, thus leading to a partition of the day made of 96 consecutive quarter-hours (U = 96 with previous notation). A typical row of the adjacency matrix we analyzed, looks like:

Conclusions

We proposed a non-stationary evolution of the latent block model (LBM) allowing us to simultaneously infer the time structure of a bipartite network and cluster the two node sets. The approach we chose consists in partitioning the entire time horizon in fixedlength time intervals to be clustered on the basis of the intensity of connections in each interval. We derived the complete ICL for such a model and maximized it numerically, by means of a greedy search, for two different networks: a simulated and a real one. The results of these two tests highlight the capacity of the model to capture non-stationary time structures.

Figure 1 :

 1 Figure 1: Aggregated number of interactions between conference attendees per quarter-hour.

  i =k and similarly for c j and yu .

  and a factorizing Dirichlet a priori distribution to Φ: ν(Φ|K, D) = Dir K (ω; α, . . . , α) × Dir D (β; γ, . . . , γ).

  conference attendees 52 and 26, between 9am and 9.15am have spoken for 16 × 20s ≈ 5m30s.The greedy search algorithm converged to a final ICL of -53217.4, corresponding to 23 clusters for nodes (people) and 3 time clusters. In Figure (2a) we show how daily quarterhours are assigned to each cluster: the class C 1 contains intervals marked by a weaker intensity of interactions (on average), whereas intervals inside C 3 are characterized by the highest intensity of interactions. This can either be verified analytically by averaging estimated Poisson intensities for each one of the three clusters or graphically by looking at Figure(2b). In this Figure we computed the total number of interactions between conference attendees for each quarter-hour and it can clearly be seen how time intervals corresponding to the higher number of interactions have been placed in cluster C 3 , those corresponding to an intermediate interaction intensity, in C 2 and so on. It is interesting to remark how the model can quite closely recover times of social gathering like the lunch break (13.00-15.00) or the "wine and cheese reception" (18.00-19.00). A complete program of the day can be found at: http://www.ht2009.org/program.php. Connections for every time interval.

Figure 2 :

 2 Figure 2: The aggregated connections for every time interval (2b) and time clusters found by our model (2a) are compared.

T and U are linked by the following relation: T = ∆ u U .

Self loops are considered here for seek of simplicity. The model can easily be extended to graphs with undirected links and/or no self loops.

Greedy search algorithms are path dependent and they could converge to local maxima.

More informations can be found at: http://www.sociopatterns.org/datasets/hypertext-2009-dynamic-contact-network/.