Thomas Amoretti
email: thomas.amoretti@ensam.eu

Florian Huet

Pierre Garambois

Lionel Roucoules

Configurable dual rotor wind turbine model based on BEM method: Co-rotating and counter-rotating comparison

Keywords: Wind-turbine, Counter-rotating, Co-rotating, Blade Element Momentum, Wake modelling, Dual-rotor Aerodynamics

de documents 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.

Introduction

Wind turbine technology is an interesting power technology in the perspective of the fight against global warming. Indeed, its emission factor is one of the lowest among all power production technologies (about 11 gCO 2eq /kW h over its life cycle against 24 gCO 2eq /kW h for hydropower, 48 gCO 2eq /kW h for solar PV and 490 gCO 2eq /kW h for gas) [START_REF] Schlömer | Technology-specific cost and performance parameters[END_REF]. Despite this advantage, wind energy production still represents a relatively small share of total electricity production (5.3 % in 2019) [START_REF] Iea | Key world energy statistics[END_REF]. In this context, research efforts should focus on improving performances of wind energy.

The performance of a wind turbine is described through its power coefficient (ratio of extracted power to wind power) curve. An example of such curve is displayed in Figure 1, where λ is the tip speed ratio (ratio between linear speed at the blade tip and upstream wind speed). To improve the performance of wind turbine technology, some solutions such as tip rotors, coplanar multiple rotors, shrouded wind turbine and unconventional power transmissions systems have been investigated [START_REF] Watson | Future emerging technologies in the wind power sector: A european perspective[END_REF]. We will focus on a specific technology, which is the coaxial multiple rotor wind turbines. Dual rotor wind turbine are an innovative technology that was developed since early 2000s [START_REF] Oprina | A review on counter-rotating wind turbines development[END_REF]. Such configuration allows recovering the outgoing flow at the back of the first turbine and use a common structure for both rotors. B. Newman [START_REF] Newman | Multiple actuator-disc theory for wind turbines[END_REF] demonstrated that the theoretical maximal power coefficient can rise from 0.594 for a single rotor (Betz limit) [START_REF] Betz | The maximum of the theoretically possible exploitation of wind by means of a wind motor[END_REF], to 0.64 according to the actuator theory applied to two rotors wind turbines.

Several studies focus on experimental performance investigation of counterrotating wind turbines [START_REF] Herzog | Performance and stability of a counter-rotating windmill using a planetary gearing: Measurements and Simulation[END_REF][START_REF] Adaramola | [END_REF][START_REF] Wang | An experimental study on the aerodynamic performances and wake characteristics of an innovative dualrotor wind turbine[END_REF]. They show that the extracted power of dual-rotor wind turbine can increase from 7% up to 12% compared to the single rotor system. Experiments have shown that using multiple rotor wind turbines in counter-rotating configuration leads to a better energy production than in co-rotating configuration [10; 11]. These experiments demonstrated that the performances of counter-rotating configuration are improved up to 20 % compare to co-rotating performances.

The fundamental question is how to define the technical parameters of both rotors of the wind turbine to approach the theoretical limit. In order to answer this question, we propose a model of dual rotor wind turbine based on BEM theory. Such model can be used as a design tool and allows optimization process for dual rotor wind turbine.

This work presents a modified BEM theory for co-rotating and counterrotating wind turbines. This latter theory takes into account the influence of the first rotor on the wake axial and tangential velocities, the velocities' evolution depending on the distance between the two rotors and the rotation direction of the second rotor. The model takes various geometrical inputs, such as the blade length, the number of blades, the chord, twist angles and profiles distribution along the blade. These parameters can vary for both rotors. Physical properties, such as wind velocities (axial and tangential), air density and dynamic viscosity are also considered. Finally, the rotational speed is also an input parameter. Technically, this latter can be controlled through pitch control or power electronic interface [START_REF] Apata | An overview of control techniques for wind turbine systems 10 e00566[END_REF]. The outputs of the model are local wind speed, axial and tangential forces along the blade and therefore power output of the whole dual rotor turbine.

The paper is organized as follows: the first part review the state of the art of dual rotor modelling. The second part address the new dual-rotor BEM model. The third part presents a reference case study. The last part compares the results of the model developed with experimental and numerical data extracted from the literature and investigates multiple rotors wind turbine performances. The results for co-rotating and counter-rotating configurations are compared.

State of the Art

Dual rotor aerodynamic models

Various methods exist to model the aerodynamic behavior of dual rotor wind turbine, including Computational Fluid Dynamics (CFD) models [START_REF] Kumar | Computational and experimental analysis of a counter-rotating wind turbine system[END_REF][START_REF] Koehuan | Investigation of counter-rotating wind turbine performance using computational fluid dynamics simulation[END_REF][START_REF] Rezaeiha | Wake interactions of two tandem floating offshore wind turbines: CFD analysis using actuator disc model[END_REF][START_REF] Faisal | Aerodynamic performance and flow structure investigation of contra-rotating wind turbines by CFD and experimental methods[END_REF] or Vortex Lattice Method (VLM) [START_REF] Lee | Velocity interference in the rear rotor of a counterrotating wind turbine[END_REF][START_REF] Rosenberg | A prescribed-wake vortex lattice method for preliminary design of co-axial, dual-rotor wind turbines[END_REF][START_REF] Jeon | Unsteady aerodynamics of offshore floating wind turbines in platform pitching motion using vortex lattice method[END_REF]. These both methods have a high calculation cost [START_REF] Wang | State of the art in the aeroelasticity of wind turbine blades[END_REF].

Another modelling technique is the Blade Element Momentum (BEM) method, which consists of dividing the rotor blades into a series of segments and analyzing the aerodynamic forces acting on each segment. The BEM theory is an efficient and a very fast modelling method used in many studies concerning single rotor wind turbine [START_REF] Wang | State of the art in the aeroelasticity of wind turbine blades[END_REF][START_REF] Lanzafame | Fluid dynamics wind turbine design: Critical analysis, optimization and application of BEM theory[END_REF][START_REF] Plaza | Comparison of BEM and CFD results for MEXICO rotor[END_REF][START_REF] Bavanish | Optimization of power coefficient on a horizontal axis wind turbine using bem theory[END_REF][START_REF] Rahimi | Investigation of the current yaw engineering models for simulation of wind turbines in BEM and comparison with CFD and experiment[END_REF][START_REF] Martínez | An improved BEM model for the power curve prediction of stall-regulated wind turbines[END_REF][START_REF] Bhargava | A case study of wind turbine loads and performance using steady-state analysis of[END_REF][START_REF] Dai | Aerodynamic loads calculation and analysis for large scale wind turbine based on combining BEM modified theory with dynamic stall model[END_REF]. As a fast computational technique, it can quickly test several configurations and pave the way for optimization. Therefore, we will focus on BEM modelling in this study.

The implementation of the BEM theory for the dual rotor configuration requires adaptations to the method. The inflow velocities for each element of the second rotor are totally modified axially and tangentially compared to the basic BEM theory. Only a few authors propose modelling approaches based on the BEM for counter-rotating wind turbine.

First, S.N Jung et al. [START_REF] Jung | Aerodynamic performance prediction of a 30kw counter-rotating wind turbine system[END_REF] investigates the performances of a 30kW counter-rotating model using a BEM model. The influence of the first rotor on the second one is taken into account from experimental data by interpolating considered axial flow velocity measured at several distances behind the first rotor. The influence of the swirl induced by the first rotor in the wake is not considered for the velocities input in the second rotor.

S. Lee et al. [START_REF] Lee | Effects of design parameters on aerodynamic performance of a counter-rotating wind turbine[END_REF] proposed a BEM model for counter rotating wind turbine. The main assumptions of this paper are that the rear rotor is operating in the fully developed stream tube of the first rotor, which leads to a reduction of the inflow wind speed by a factor of (1-2a) for the rear rotor. Here a is the axial induction coefficient for a specific location on the first rotor's blade. The tangential velocity induced by the front rotor in the inflow of the second one is not taken into account.

B. Hwang et al. [START_REF] Hwang | Optimization of a counter-rotating wind turbine using the blade element and momentum theory[END_REF] published another study based on the same principle but with a distance between the two rotors fixed at 0.33 D, where D is the diameter of the rotor. In order to consider the fact that the second rotor is not in the fully developed stream tube, the speed reduction factor is set to (1-1.6a) according to the study [START_REF] Lee | Velocity interference in the rear rotor of a counterrotating wind turbine[END_REF]. Like previous studies, there is no use of the tangential velocity induced by the front rotor.

F. Yin et al. [START_REF] Yin | A blade element momentum model for dual-rotor wind turbines considering inter-rotor velocity interferences 2265[END_REF] published a study of a BEM theory that calculates both interactions between the two rotors for the axial induction part. For the tangential component, only the action of the first rotor on the second one is considered. The impact of the distance between the turbines on the airflow velocities is not taken into account. Finally, F. Mühle et al. [START_REF] Mühle | The effect of rotational direction on the wake of a wind turbine rotor -a comparison study of aligned co-and counter rotating turbine[END_REF] published an article comparing a wind tunnel test of co and counter-rotating wind turbines performances with calculated performances from a BEM theory applied to this configuration. They take into account the effects between the two rotors by integrating wind speeds measured data as entries for the BEM model.

Discussion and objectives

After reviewing the literature (see table 1), modelling of dual rotor wind turbine with CFD or VLM methods are time expansive. BEM method is faster, but most studies do not take into account tangential effects of the flow that are essential to separate co-rotating and counter-rotating performances. Besides, there is no consideration of the distance between the two rotors in the literature, which is an important design parameter to configure. The present study proposes a fast, accurate and fully configurable model that can provide the mechanical power output for any distance between the two rotors and rotors configurations. [START_REF] Wang | State of the art in the aeroelasticity of wind turbine blades[END_REF] and specific BEM characteristics

Model

Dual-rotors Blade Element Momentum model

Model assumptions

The assumptions of the model are as follows :

• The flow is considered inviscid and incompressible;

• The system is in a stationary state;

• The retroaction from the second rotor to the first one is neglected;

• There is no radial speed taken into account;

• The expansion of the flow at the rear of the first rotor is not considered;

The last assumption is made because the expansion of two adjacent concentric annular sections would result in the intersection of the two respective flows and give complex flow velocity calculations.

The dual-rotor wind turbine will be modelled with two distinct rotors whose performances are calculated using the BEM theory.

Blade Element Momentum theory : single stage

The Blade Element Momentum theory is a common calculation technique applied to wind turbine. This theory is described in several books on wind turbine aerodynamics [START_REF] Gasch | Wind power plants: Fundamentals, design, construction and operation[END_REF][START_REF] Manwell | [END_REF][START_REF] Hansen | Aerodynamics of wind turbines[END_REF]. It combines the theories of momentum conservation and the blade element aerodynamic forces to calculate the power extracted from the flow by the rotor. The figure 2 represents a single wind turbine rotor with a radius R prop with a rotational speed of Ω. R hub represents the hub radius. The radial, tangential and axial directions are identified by vectors (⃗ e r , ⃗ e θ , ⃗ e x ).

Power calculation

The blade is discretized into a finite number of N blade elements. Local forces and torques are calculated on each element, according to the aerodynamics effects. The input flow is divided in concentric annular section that correspond to each element. An element is represented at a distance r from the rotational axis ( ⃗ e x ). The geometry of the element is described in figure 3 where c is the chord length of the element, V ax0 is the flow velocity far upstream the turbine (fully axial) and dr is the length of the element. The element is also defined by two aerodynamic profiles that can be different at element's foot and element's tip. The flow velocities and forces applied on the blade element located at a distance r from the center of the turbine are detailed in the Figure 4. The angle between the chord axis and the total speed relative to the blade (V br ) is the attack angle, α. The angle between the rotation plane and the chord axis is the twist angle β. The angle Φ is the sum of the two previous angles. The aerodynamic forces applied to the profile are the drag force ( ⃗ F D ) in the direction of ⃗ V br and the lift force ( ⃗ F L ) in the perpendicular direction at ⃗ V br . ⃗ F θ and ⃗ F x are the projections of the previous forces respectively on ⃗ e θ and ⃗ e x axis.

The forces applied on an element are calculated with equations 1 and 2 where ρ is the air density.

dF x = 1 2 .ρ.c.dr.V 2 br .[C L .cos(Φ) + C D .sin(Φ)] (1)
dF θ = 1 2 .ρ.c.dr.V 2 br .[C L .sin(Φ) -C D .cos(Φ)] (2)
In equations 1 and 2 C L and C D are aerodynamic coefficients (respectively lift and drag coefficients). We determine these later coefficients according to the attack angle (α), the local Reynolds number and the airfoil shapes composing the blade element. Good aerodynamic data concerning these coefficients is an essential prerequisite for the BEM calculation.

The torque developed by one element is described in equation 3.

dM = r.dF θ (3)
Then the power extracted by the whole first rotor is calculated with equation 4.

P turbine = B. Rprop R hub Ω.dM (4)
With B number of blades. The power of the wind passing through the area swept by the blades (P wind ) is expressed in equation 5.

P wind = 1 2 .ρ.π.(R 2 prop -R 2 hub ).V 3 ax0 (5)
The power coefficient is defined as the ratio of the power extracted by the turbine to the total power of the wind (equation 6).

C p = P turbine P wind (6)

Wind speeds around the blade

The blade element modifies axial and tangential speeds in the surrounding flow. This effect is equivalent to new speeds induced by the blade element that add up with flow velocities. In order to calculate the forces applied on an element, we need to determine the value of the induction.

V bax is the total axial flow seen by the blade (equation 7).

V bax = (1 -a).V ax0 (7)
V bax is the superposition of the upstream flow velocity V ax0 and the axial induction effect -a.V ax0 (equation 7).

V bt is the total tangential flow seen by the blade (equation 8).

V bt = V btr + V btind = (1 + a ′ ).Ω.r (8)
V bt is the superposition of the tangential relative wind speed seen by the moving blade (V btr = Ω.r) and the effect of induction (V btind = a ′ .Ω.r) (equation 8). a and a ′ are respectively the axial and tangential induction coefficients for the first rotor. Since the tangential induction is directed in the opposite direction of the blade rotation direction, these two speeds add up.

The global resultant velocity V br is the vectorial composition of axial and tangential velocities seen (equation 9).

V br = V 2 bax + V 2 bt (9)

Induction coefficients calculation : fixed point algorithm

Equations of the Blade Element Momentum leads to the calculation of the induction coefficients, which allows on the one hand to calculate the forces applied on the profile and in the other hand to calculate the velocities in the flow just behind the wind turbine. An important issue of the BEM theory is the calculation of these induction coefficients (a and a'). We solved it through an iterative fixed-point algorithm where values of induction depends on previous induction setting until the convergence criterion (ε) is insured. This latter is described in algorithm 1.

For each n th iteration, the induction coefficients are given by equations 10 and 11.

a (n+1) = 1 4.sin(Φ) 2 s.[C L .cos(Φ)+C D .sin(Φ)] + 1 (10) a ′ (n+1) = 1 4.sin(Φ).cos(Φ) s.[C L .sin(Φ)-C D .cos(Φ)] -1 (11)
Algorithm 1 Fixed point algorithm description

Require: a (0) ← 0 and a ′ (0) ← 0 and a (1) ← (eq.10) and a ′ (1) ← (eq.11) and load aerodynamic tables 10) 11) end while

while |a (n) -a (n-1) | > ε or |a ′ (n) -a ′ (n-1) | > ε do Φ ← arctan( (1-a (n) ).V ax0 (1+a ′ (n) ).Ω.r ) α ← Φ -β C L ← aerodynamic table lift coeff(α) and C D ← aerodynamic table drag coeff(α) a (n-1) ← a (n) a ′ (n-1) ← a ′ (n) a (n) ← a (n+1) (Equation
a ′ (n) ← a ′ (n+1) (Equation
Where C L and C D are the lift and drag coefficients. s is the local solidity, defined as the ratio of the sum of blade chords at this radius over the perimeter drawn by this radius (equation 12).

s = c.B 2π.r ( 12
)
Where B is the number of blades of the rotor.

Once the induction coefficients are calculated and the relative speed seen by the element is known, it is possible to calculate the forces applied on this latter.

BEM improvements

The BEM theory is enhanced with some improvements. As it stands, this theory does not match with experimental data for high values of axial induction coefficient. To address this issue, Spera's correction is implemented [START_REF] Hansen | Aerodynamics of wind turbines[END_REF]. This latter proposes to calculate the (n + 1) th axial coefficient with equation 13 if the induction factor is above the critical value (a c ) of 0.2.

a (n+1) = 1 2 .[2 + K.(1 -2.a c ) -(K.(1 -2.a c ) + 2) 2 + 4.(K.a 2 c -1)] (13) With K = 4.sin(Φ) 2 s.[C L .cos(Φ) + C D .sin(Φ)] (14)
This equation replaces equation 10 in the case of high induction.

Another correction takes into account the finite number of blade and the tip losses [START_REF] Hansen | Aerodynamics of wind turbines[END_REF]. This introduces a so-called Prandtl's tip loss factor that modifies the induction coefficients for elements close to the tip of the blade (see equations 15 and 16). The closer the elements are to the blade tip, the more their performances are degraded.

F = 2 π × arccos(e -f ) (15)
With

f = B 2 × R prop -r r.sin(Φ) (16)
Prandtl's tip loss factor modifies respectively equations 10 and 11 by equations 17 and 18.

a (n+1) = 1 4.F.sin(Φ) 2 s.[C L .cos(Φ)+C D .sin(Φ)] + 1 (17) a ′ (n+1) = 1 4.F.sin(Φ).cos(Φ) s.[C L .sin(Φ)-C D .cos(Φ)] -1 (18)
In case of high inductions factors, equation 14 becomes 19.

K = 4.F.sin(Φ) 2 s.[C L .cos(Φ) + C D .sin(Φ)] (19)
3.3. Calculation of the wind speeds in the wake of the rotor We need the induction values calculated during the BEM to calculate the speed modification in the wind and therefore the forces applied on the second turbine. The values of the airflow velocities just behind the rotor are respectively (1 -a).V ax0 for the axial one and -a ′ .Ω.r for the tangential one, oriented in the opposite direction of the turbine rotation direction. Far downstream, momentum theory states that these values reach respectively (1-2a).V ax0 and -2a ′ .Ω.r. This means there is an influence of the distance to the turbine for the flow velocities. The flow velocities evolution is displayed in Figure 5a for the axial speed and Figure 5b for the tangential speed. The turbine is in the 0 position on the x -axis. These speed evolutions do not We can define the velocity evolution over the x -axis with a distance coefficient (C distance ). The distance coefficient should have the value 0 on x → -∞, 1 when x = 0 and 2 when x → +∞. The evolution around x = 0 should be linear [START_REF] Hansen | Aerodynamics of wind turbines[END_REF]. In this work, we use the relation proposed by O. Gur [START_REF] Gur | Extending blade-element model to contra-rotating configuration 638[END_REF] for propeller stages (equation 20). where x is the distance from the turbine, R prop is the turbine radius and θ the angle formed by the x-axis and the straight line passing through the blade tip (see figure 6).

C distance (x) = 1 + cos(θ) = 1 + x x 2 + R 2 prop (20) 0 x Rprop θ ex
The values of velocities in the flow can be described by equations 21 and 22. These velocities are calculated for each concentric annular element, so also have a radial dependence. With this notation V ax0 is V ax (-∞).

V ax (x) = V ax0 .(1 -C distance (x).a) (21) V t (x) = -a ′ .Ω.r.C distance (x) (22)

Application of the BEM on a second rotor

In the case of dual-rotor wind turbine, the BEM calculation is carried out in the same way on the second rotor. The only changes are the inflow modification induced by the presence of the first rotor and the distance between the two turbines. The flow modification has to be considered in the calculation of the performances of the second rotor. That supposes a new tangential entry that does not exist in the first rotor case. This new rotational input is generated by the first rotor, the latter also generating an axial deceleration. Figure 7 represents the whole dual-rotor system. Here, the index 1 for each variable refers to the first rotor when the index 2 stands for the second rotor. For the second rotor, located at a distance x = d from the first one, the procedure of Blade Element Momentum calculation remains mainly the same but with the following changes :

• First, the wind upstream velocities considered are now V ax (d) and V t (d).

The expression of global axial flow seen by the blade is described in equation 23.

V bax2 = (1 -a 2 ).V ax (d) (23)
• Second, the global flow seen by the blade in the tangential direction depends on the rotation direction of the blade. The tangential relative wind speed is described in the equation 24.

V bt2 = V t (d) + V btr2 + V btind2 (24)
The global tangential speed without the induction effect is the superposition of the wind seen by the rotation of the blade (V btr2 = Ω 2 .r) and the upwind tangential wind velocity at the d location on the x-axis (V t (d)). Figure 8 represents both co-rotating and counter-rotating configurations. If the blade of the second rotor turns in the opposite direction compare to the first one (counter-rotating case, see figure 8a), the tangential wind inflow (V t (d)) is in the same direction as the movement of the blade. The tangential wind inflow is therefore subtracted from the relative tangential flow seen by the rotating blade (V btr2 ). In the co-rotating case (see figure 8b), the tangential inflow is in the opposite direction of the movement of the second blade, therefore the tangential relative wind speed and the tangential inflow speed add-up.

The tangential induction effect applies in the whole tangential speed seen by the blade. By considering the direction of rotation, we have the total tangential speed for the counter-rotating configuration in equation 25. while for the co-rotating configuration, the total tangential speed is described in equation 26.

V bt2 (counter) = (1 + a ′ 2 ).(Ω 2 .r -V t (d)) (25)
V bt2 (co) = (1 + a ′ 2 ).(Ω 2 .r + V t (d)) (26)
The other equations of the BEM remains unchanged, especially the equations to calculate the induction coefficients. Only the inflow wind speed vector (norm and angle) is modified. We apply on the second rotor the same procedure that the one described in the BEM part. The final power of the system and power coefficient is the sum of powers extracted by the two rotors (equation 27).

P turbinetot = P turbine1 + P turbine2 (27)
The final power coefficient (equation 28) is the ratio between the total power and the total power of the wind (see equation 5).

C ptot = P turbinetot P wind (28)

Reference model

The global model of dual-rotor wind turbine system was presented in the previous section. This part present our object of study to apply this model.

Global settings

The large number of modifiable parameters involved in the BEM calculation and the geometry description of the rotor (radius, twist angles, airfoils...) requires a precise definition of the wind turbine used to allow reproducibility of the tests. F. Mühle et al. [START_REF] Mühle | The effect of rotational direction on the wake of a wind turbine rotor -a comparison study of aligned co-and counter rotating turbine[END_REF] proposes an experimental and BEM model study of dual-rotor wind turbine with variable distance between the two rotors. Their model takes experimental input for the second rotor. We will therefore take Mühle's study as a reference to compare our fully parametric model.

In F. Mühle et al. study, both turbines have the same geometrical configuration (e.g. Rprop 1 = Rprop 2 ). Only the direction of rotation can change for the second turbine. Table 2 summarizes items that define the wind turbine used.

A wind turbine of this size (R prop = 0.4425m) is explained by the fact that experimental wind tunnel tests limit the size of prototypes. This corresponds to a domestic wind turbine with an output of a few hundred watts. The blade used in this study come from the work of P. Krogstad and J. Lund [START_REF] Krogstad | An experimental and numerical study of the performance of a model turbine[END_REF] where the blade specifications are fully described. Geometrical specifications of the blade are summarized in the table 4.

We used the same geometrical specifications and wind conditions in our model.

Airfoil description

The determination of lift and drag coefficients depending on the attack angle and the local Reynolds number is essential to calculate the induction coefficient and apply the BEM theory (see Algorithm 1). Figure 9 describes aerodynamic performances of the S826 profile depending on the attack angle for different Reynolds' numbers. The aerodynamic data to describe the S826 airfoil performances was generated with Qblade software. Lift and drag coefficient calculation for attack angles from -20°to 20°are made with a transition amplification ratio of N crit = 3 which correspond to a dirty wind tunnel condition [START_REF] Drela | Xfoil user guide[END_REF]. This ratio set the location of the laminar-turbulent transition in the boundary layer, which has an impact on aerodynamic performances. The curve in then extrapolated to 360°with Viterna's method [START_REF] Mahmuddin | Airfoil lift and drag extrapolation with viterna and montgomerie methods[END_REF].

We set the number of blade elements to 100 for each blade to ensure convergence of the BEM results. The convergence criterion of the fixed point algorithm is set to ε = 10 -3 .

Calculation protocol

The power coefficients calculations were made with the same protocol as described in the work of F. Mühle et al. [START_REF] Mühle | The effect of rotational direction on the wake of a wind turbine rotor -a comparison study of aligned co-and counter rotating turbine[END_REF]. The rotational speed of the first wind turbine is kept constant to have a tip speed ratio (λ 1 ) of 6. The tip speed ratio is defined as the ratio of the linear speed at the blade tip to the upstream axial wind speed (equation 29). Considering the blade length and the wind upstream velocity V ax0 , the first rotor rotates at Ω 1 = 1300 rpm.

λ 1 = Ω 1 .R prop1 V ax0 = 6 (29)
In our application case, only the second rotor's tip speed ratio variate with the change of the rotational speed of the second rotor (Ω 2 ). The tip speed ratio of the second rotor (λ 2 ) is defined by keeping the same upstream axial wind speed (equation 30).

λ 2 = Ω 2 .R prop2 V ax0 (30)
Three distances between the two rotors were taken into account. These distances are 4 × R prop , 7 × R prop and 10.3 × R prop .

For every distance and every tip speed ratio, the performances of both counter-rotating and co-rotating configuration were calculated. The global performance of a configuration is the sum of the two power coefficients of the two rotors (see equations 27 and 28).

Results and discussion

In this section, we use the case study described in section 4 to validate our model. We will exploit our model by comparing counter-rotating and corotating configurations and single and double rotors cases. We will compare to F. Mühle et al. [START_REF] Mühle | The effect of rotational direction on the wake of a wind turbine rotor -a comparison study of aligned co-and counter rotating turbine[END_REF] study, mentioned as "reference". The present study will be referred as "current study".

Accuracy of the model 5.1.1. Validation of BEM on a single stage wind turbine

We compare the results of our BEM model on a single stage wind turbine with the experimental and BEM results of the reference. The comparison is carried out for λ = 1 to λ = 7 as this is the range of λ used for the dual rotor investigation in the reference.

Figure 10 represents the C p calculated by the present study compare to experimental and BEM calculation of the reference.

Table 5 shows the relative differences between current BEM study and experimental results and between reference BEM study and experimental results.

The comparison shows that the average difference between the BEM model and the experiment is comparable between the current study and the reference (around 19% error in average). The standard deviation of relative error is higher for current study, but high C p values' region (from λ = 4 to λ = 7) is described with higher accuracy.

Dual rotors models comparison results

The results for the calculations associated to the case study are depicted in figures 11, 12 and 13 for the three different distances between the two rotors and compared to the results of the reference [START_REF] Mühle | The effect of rotational direction on the wake of a wind turbine rotor -a comparison study of aligned co-and counter rotating turbine[END_REF]. Every figure depicts the performances of co-rotating and counter-rotating configurations against λ 2 for the experimental reference, the model reference and the model of the present study.

Figure 11 displays results for a 4 × R prop distance between the two rotors. Figure 12 presents results for a 7 × R prop distance between the two rotors.

Figure 13 provides results for a 10.3 × R prop distance between the two rotors. The relative difference between BEM studies and experimental results was calculated for each point.

The means, the standard deviations and the extreme values of all the above-mentioned relative difference points series are summarized in Table 6.

Table 6 shows that the model developed in this study is in good accordance with experimental data for low distances between rotors, with a mean relative error of 5.1 % for d = 4 × R prop in counter-rotating case and 5.3 % in co-rotating case. For longer distances, the relative error of the model increases to 15.3 % for d = 10.3 × R prop concerning the counter-rotating case and 16.4 % for co-rotating case. This error is reasonable as the model is not fed by experimental data for second rotor inputs contrary to the reference BEM. The accuracy of calculated BEM of the reference appears to stay constant with the growing distance in both configurations.

Discussion about models comparison

In the subsection 5.1.2, we can see that the smaller the distance between the two turbines, the more accurate the BEM calculation. We will now discuss the reason that explain the difference between the curves. The viscosity of the flow outside the current tube crossing the wind turbine system is neglected in our BEM calculation. The flow outside the current tube enclosed by the rotors has an axial speed of V ax0 which is higher than the flow downstream of the first rotor. The viscosity of the fluid induces an energy transfer from the outside to the inside of the current tube, which means acceleration of the flow between the two rotors. This effect makes the flow passing the turbine to reach the velocity of V ax0 far behind the turbine in real condition. This effect is not considered in the distance coefficient based on momentum conservation theory (see figure 5). This is obviously taken into account in the BEM proposed by the reference because of the feeding of their model by wind tunnel experimental data. The further away the second rotor is from the first turbine, the greater the underestimation of the wind speed observed by the second turbine. This effect also occurs on the tangential speed. In real conditions, zero tangential velocity is expected far downstream of the rotor. The greater the distance, the greater the overestimation of the tangential wind speed observed by the second turbine. These effects tend to underestimate the C p value compare to the experiments with a bigger effect with the growing distance. In the aim of having counter-rotating wind turbine, the distance between the two rotor is small for technical and mechanical purposes. We can therefore assume that the model is accurate for our case study of counter-rotating wind turbine. F. Mühle [START_REF] Mühle | The effect of rotational direction on the wake of a wind turbine rotor -a comparison study of aligned co-and counter rotating turbine[END_REF] explains the difference between his BEM calculation and the experimental results because of the wind tunnel test conditions. It prevents the wake to expand freely and his BEM underestimate the C p calculated. On the contrary, the current study does not take into account wake expansion and is expected to overestimate the C p value compare to the model developed in the reference. These conditions are closer to the experimental setting, which takes place in a wind tunnel.

The current model shows a good agreement with experimental values in the perspective of a fast calculation and configurable model.

Co-rotating vs counter-rotating configuration comparison 5.2.1. Results of the two configurations

The figure 14 displays the performances of counter-rotating and co-rotating configurations in function of λ 2 for a distance d = 4 × R prop between the two turbines. The curves are for experimental data and BEM calculation. In our BEM simulation, for values of λ 2 < 2 we can see that the C p value is slightly higher for co-rotating configuration (about 0.5%). For λ 2 > 2 the counter-rotating configuration become clearly more efficient than co-rotating (up to 4.6 %). In the experimental results, the same behavior occurs with a switch for λ 2 = 1.5. The counter-rotating configuration appears to have a C p value 4.3 % higher than co-rotating.

Table 7 displays mean values of the global relative wind speed seen by the second blade (V br2 ) and the mean value of attacks angles (α). This table is for both rotor configuration for λ 2 = 1 and λ 2 = 3.5.

Table 7 shows a higher average value for relative wind speed in co-rotating case for both tip-speed values compare to counter-rotating value. There is a lower attack angle for co-rotating configuration compare to counter-rotating configuration for both λ 2 values. The C p value is approximately the same for both configurations for λ 2 = 1 but is higher for co-rotating case for λ 2 = 3.5. From these results we can see that in this configuration (same blade front and back rotor) the counter-rotating configuration is better in most λ 2 cases and co-rotative configuration is slightly better for low values of λ 2 .

Counter-rotating Co-rotating

λ 2 = 1 mean V br2 [m/

Discussion on the configurations

The rotational direction of the second turbine has the effect of adding or subtracting the inflow tangential component for each blade element of this turbine. This modification of the tangential speed seen by the blade has two effects :

• The first one is the modification of the resultant wind speed norm V br2 (equation 9);

• The second effect is the modification of the attack angle seen (α);

Both of these effects have an impact on the forces applied on the blade, and particularly on the tangential force useful for the power coefficient calculation (see equations 1 and 2 and figure 8). In these equations C L and C D coefficients and Φ angle depend on the attack angle and the force is proportional to V 2 br2 . The Table 7 shows that, in average, the resultant wind speed is higher for a co-rotating wind turbine, but the attack angle is lower than for the counter-rotating configuration. In the calculation of tangential force for counter-rotating case, there is a conflict between the diminution of V br2 and the evolution of the (sin(Φ).C L -cos(Φ).C D ) coefficient with the increase of α value compared to the co-rotating case. For λ 2 = 3.5 this latter coefficient is more beneficial for the counter-rotating configuration than the reduction of V br2 compare to co-rotating case as C p value for counter-rotating is higher than C p co-rotating. The higher value of the (sin(ϕ).C L -cos(ϕ).C D ) coefficient for counter-rotating case can be explained as the mean attack angle is closer to the best lift-to-drag ratio attack angle that is around α = 6 • (see figure 9).

In the case λ 2 = 1, the effects are the same except a higher α leads to a decrease of lift-to-drag ratio, the maximum being exceeded. The average attack angle is closer to the optimal angle for maximum lift-to-drag ratio in the co-rotating configuration in this case. Coupled with the increase of V br2 value, it explains why co-rotating turbines have better performances for λ 2 < 2.

Co-rotating configuration has a higher resultant wind speed norm (V br2 ) value compare to counter-rotating configuration. With appropriate attack angle, it has the potential to exceed counter-rotating performances.

Dual-rotor vs single rotor comparison 5.3.1. Results of both systems

We want to compare the performances of single and dual rotor configurations. As the two rotors can rotate independently, the tip speed ratio λ is not defined for the system. We therefore calculate C p values for different upstream wind speeds and different rotational speeds. We base ourselves on an upstream wind vector V ax0 ranging from 1 to 25 m/s with a wind step of 1 m/s. We take the following assumptions :

• The rotational speed of the single wind turbine can change between two different upstream wind speed.

• The rotational speeds of front and back rotor in a dual rotor case can vary independently of each other for every upstream wind speed.

For each wind speed, we have to determine the best functioning point for this particular wind speed, i.e. the best vector ([V ax0 ,Ω 1 ] for single rotor and [V ax0 ,Ω 1 ,Ω 2 ] for dual rotor) that maximizes the C p value. We set a rotational speed limit as Ω max = 2000 rpm to ensure the flow seen by the blade is incompressible. Unlike the previous paragraphs, where only the rotational speed of the second rotor can vary, we have in this paragraph two (or three) degrees of freedom.

For the single rotor case, we scan different rotation speed and take the best C p value. For instance, figure 15 depicts the C p curves in function of Ω 1 for five V ax0 values going from 1 m/s to 25 m/s.

In figure 15 we can see that the maximum C p value is not reached at the same Ω value for every V ax0 . The maximum C p value is not constant with the upstream wind speed (C pmax ranges from 0.18 up to 0.50). The C p can reach negatives values in the case of low V ax0 . That means that the required speed is too high in relation to the energy contained in the wind flow. For the dual-rotor configurations, we scan all possible pair values for Ω 1 and Ω 2 for every upstream wind speeds. For instance, figure 16 represents the surfaces obtained for wind speed of V ax0 = 10 m/s and V ax0 = 15 m/s. Each wind speed is represented for respectively counter-rotating (figures 16a and 16b) and co-rotating configurations (figures 16c and 16d). The distance between the two turbines is d = 4 × R prop .

In figure 16 we can see that the C p value reaches a maximum for different (Ω 1 ,Ω 2 ) pairs and the maximum values are different. For instance, in counterrotating case, for V ax0 = 10 m/s we have C pmax = 0.5164 reached for (Ω 1 ,Ω 2 ) = (1192,505) rpm. For V ax0 = 15 m/s we have C pmax = 0.5315 reached for (Ω 1 ,Ω 2 ) = (1798,808) rpm.

The performance curves for all wind speeds are formed by the best C p for each upstream wind speed for each configuration. Figure 17 shows the C p curves depending on wind speed V ax0 for single rotor and dual-rotors configurations. Examples of C p points on the curve are shown (V ax0 = 5 and 10 m/s for single rotor and V ax0 = 10 and 15 m/s for counter-rotating configuration).

In figure 17 we can see that for every wind speed, dual-rotor configurations are better than single rotor one. The max C p value for single rotor is 0.500 reached for V ax0 = 16 m/s. The max C p for counter-rotating configuration is 0.532 reached for V ax0 = 16 m/s. We can see three main areas for all configurations with an increase in C p values to V ax0 = 5 m/s, a plateau area

Discussion on performances

As shown in figure 17, the counter-rotating configuration has better performances than co-rotating for every wind speed. Which means that independently of the location of the wind turbine and thus the wind distribution, the counter-rotating configuration will provide the most energy. The decrease of the C p curves can be explained by the limitation of the rotational speed to keep the flow seen by the blade incompressible. It is impossible for the blade For some upstream wind speed, there is a difference in the optimal rotational speed between the single wind turbine and the rotational speed of the first rotor in the dual-rotor case. For example, for V ax0 = 11 m/s, the optimal rotational speed is 1354 rpm for the single rotor and 1232 rpm for the counter-rotating case (C p = 0.517). In this example, the maximum C p value for a Ω 1 set at 1354 rpm would have been C p 0.516. As a reminder, the second rotor does not have retroactive effect on the first rotor. This means that for some upstream wind speed, it is more beneficial for the whole system to depreciate the first rotor's performance to allow a bigger gain on the second one.

Conclusion and perspectives

A BEM model for dual-rotor wind turbine has been developed. It was compared to experimental data and to a reference BEM model fed by experience for several tip speed ratios (of the second rotor λ 2 ) and distances between the two rotors. The results show a good performance match for co-rotating and counter-rotating configurations.

Comparison between co and counter rotating configurations shows that counter-rotating can produce a C p value 4.3% higher in case of same front and back rotors. Both dual rotor configurations produce higher C p value than single rotor for different wind speeds investigated (10.6 % in average for counter-rotating configuration).

The model developed is fully configurable from wind turbine geometrical aspects (size of blades, number of blades, size of the hub) to local blade element setting distribution (twist angle, chord length, profile). This setting can be different for both upstream and downstream rotors. Wind flow conditions are also configurable (wind speed, air density, dynamic viscosity). The second rotor can be either in co-rotating or counter-rotating configuration. The effect of distance between the two rotors is taken into account. The model does not require experimental data to provides results on power extracted by the rotors.

The computational time with an 11th Gen Intel © Core™ i7-1165G7 @ 2.80GHz/1.69 GHz and 16 Go RAM, using Matlab software is 0.3 s for dual rotor calculation of one single C p value. A great advantage of a such fully configurable, accurate and fast calculation model is to explore a wide design space. This enables the application of optimization algorithms in order to improve the performances of wind turbines.

Acknowledgment

The French region SUD (Provence-Alpes-Côte-d'Azur) supported this work. The authors gratefully acknowledge this support.

Figure 1 :

 1 Figure 1: Performance curve of a wind turbine

Figure 2 :

 2 Figure 2: Wind turbine parameters definition

Figure 3 :

 3 Figure 3: Blade element geometry

Figure 4 :

 4 Figure 4: Flow velocities and forces seen by an element

Figure 5 :

 5 Figure 5: Axial (a) and tangential (b) wind speed evolution along the x-axis

Figure 6 :

 6 Figure 6: Distance coefficient angle

Figure 7 :

 7 Figure 7: Dual rotors wind turbine definition

Figure 8 :

 8 Figure 8: Second rotor's velocities in counter rotating (a) and co-rotating (b) configuration

Figure 9 :

 9 Figure 9: S826 profile aerodynamic performances with Ncrit = 3

Figure 10 :

 10 Figure 10: C p calculation comparison between present study and reference study : single stage wind turbine

Figure 11 :

 11 Figure 11: Co and counter-rotating performances for a distance of 4Rprop

Figure 12 :

 12 Figure 12: Co and counter-rotating performances for a distance of 7Rprop

Figure 13 :

 13 Figure 13: Co and counter-rotating performances for a distance of 10.3Rprop

Figure 14 :

 14 Figure 14: Comparison of counter-rotating and co-rotating performances

Figure 15 :

 15 Figure 15: Single rotor performances depending on the wind speed and the rotational speed

Figure 16 :

 16 Figure 16: C p surfaces for counter-rotating and co-rotating configurations

Figure 17 :

 17 Figure 17: Performances depending on the wind speed

Table 1 :

 1 Aerodynamic models comparison from

		Speed Accuracy	BEM specific criteria
				tangential speed distance
	CFD based	+	++++
	VLM based	+++	+++
	BEM based (literature) ++++	+++	rare	no
	Present study	++++	+++	yes	yes

Table 3

 3

	summarizes items that define the wind tunnel test conditions.

Table 2 :

 2 Characteristics of the wind turbine case study

	Characteristic	Value	Symbol	Unit
	Upwind wind speed	10	V ax0	[m/s]
	Air density	1.225	ρ	[kg/m 3 ]
	Air dynamic viscosity 1.809.10 -5	µ	[P a.s]

Table 3

 3

		: Wind conditions
	radius position [m]	0.068 0.082 0.098 0.113 0.127	0.143 0.158
	chord (c) [m]	0.081 0.080 0.077 0.073 0.069	0.065 0.061
	twist angle (β) [rad] 0.647 0.568 0.501 0.441 0.391	0.349 0.315
	radius position [m]	0.172 0.188 0.203 0.217 0.233	0.248 0.262
	chord (c) [m]	0.058 0.054 0.051 0.048 0.046	0.044 0.042
	twist angle (β) [rad] 0.285 0.256 0.228 0.206 0.188	0.171 0.155
	radius position [m]	0.278 0.293 0.307 0.323 0.338	0.352 0.368
	chord (c) [m]	0.040 0.038 0.036 0.035 0.033	0.032 0.031
	twist angle (β) [rad] 0.139 0.127 0.115 0.103 0.093	0.082 0.072
	radius position [m]	0.383 0.397 0.413 0.428 0.442
	chord (c) [m]	0.030 0.029 0.028 0.027 0.026
	twist angle (β) [rad] 0.062 0.051 0.039 0.019 -0.013

Table 4 :

 4 Blade definition

Table 5 :

 5 Relative differences between BEM studies and experiment (single rotor) (in %)

Table 6 :

 6 Relative differences between BEM studies and experiment (dual rotor) (in %)

				Current BEM study
		counter-rotating		co-rotating
	Distance	mean std min max mean std min max
	4 × R prop	5.1	2.2 2.1	8	5.3	1.5 2.3	7.1
	7 × R prop	9.2	4.2 1.8 14.5	8.4	2.8 3.5 11.9
	10.3 × R prop	15.3 8.2 1.9 24.2 16.4 8.1 2.8 27.6
				Reference BEM study
		counter-rotating		co-rotating
	Distance	mean std min max mean std min max
	4 × R prop	10.6 3.4	4	15.1 10.3 2.6 6.2 14.5
	7 × R prop	11.2 3.3 4.3 15.2	9.6	3.1 3.9	14
	10.3 × R prop	11.2 3.1 4.9 15.4 10.6 3.3 3.4 14.9

Table 7 :

 7 Rotational direction effects comparison for different tip speed ratios

		s]	7.3	8.5
		mean α [deg]	24	16
		C p	0.488	0.489
	λ 2 = 3.5	mean V br2 [m/s] mean α [deg]	20.6 1	22.1 -1
		C p	0.498	0.491