N Lahellec

J C Michel

H Moulinec

P Suquet
email: suquet@lma.cnrs-mrs.fr

ANALYSIS OF INHOMOGENEOUS MATERIALS AT LARGE STRAINS USING FAST FOURIER TRANSFORMS

Introduction

Every mat erial is nat urally inhomogeneous at a small enough scale. T he aim of hom ogenizat ion is t o replace a microscopically inhomogeneous mat erial by a macroscopically homogeneous medium, t he propert ies of which are t he effect ive propert ies of t he init ial inhomogeneous mat erial. However, nonlinear phenomenas such as fat igue or rupt ure are oft en t riggered by local st ress or st rain concent rat ions at small scale. T herefore, besides effect ive propert ies, it is also essent ial t o invest igat e how local fields are infl uenced by, and infl uence in t urn, t he microst ruct ure of a mat erial.

P recise det erminat ion or accurat e est imat ion of local fields may involve considerable comput at ional effort s. Direct numerical simulat ions of composit es wit h periodic microst ruct ures have been performed by t he Finit e Element Met hod (FEM) since t he early sevent ies, but despit e numerous achievement s, t he applicat ion of t his met hod t o t hree-dimensional problems is st ill limit ed by t he complexit y of meshes. T he need for invest igat ing mat erials wit h complex microst ruct ure wit hout meshing provided a first mot ivat ion for developing a new comput at ional met hod, specifically devised for microst ruct ures, which could make direct use of images of microst ruct ures, as delivered by Scanning Elect ron Microscopy or Comput ed Tomography ( [START_REF] Moulinec | A fast numerical met hod for comput ing t he linear and nonlinear propert ies of composit es[END_REF] [2] [START_REF] Michel | Effect ive propert ies of composit e mat erials wit h periodic microst ruct ure: a comput at ional approach[END_REF]).

A second mot ivat ion for developing a met hod capable of simulat ing easily t he response of different microst ruct ures st ems from t he recent t heoret ical models of t he effect ive propert ies of nonlinear composit es (see t he review art icle P ont e Cast a ñeda and Suquet [START_REF] Castañeda | Nonlinear composit es[END_REF] and t he references herein). Direct comparison of t hese models wit h experiment al dat a is oft en difficult , in t hat many different phenomena may occur simult aneously in t he nonlinear regime (plast icit y but also damage, int erphase debonding, grain boundary sliding...). A comparison wit h numerical simulat ions, alt hough only a st ep t owards a complet e validat ion, permit s t o " act ivat e"all nonlinear mechanisms independent ly.

2.

Solving a periodic elasticity problem with the FFT As a reminder for readers who are not familiar wit h t he FFT met hod, t he case of linear elast icit y at small st rains is discussed first . T he microst ruct ure of t he composit e mat erial is described by a represent at ive volum e elem ent (r.v.e.) V comprised of N homogeneous phases. T he N individual const it uent s are linear elast ic, wit h st iffness c(x) and perfect ly bonded accross t heir int erfaces. T he volume element is subject ed t o an average st rain E and t he local problem t o be solved for t he local st ress and st rain fields reads as

σ(x) = c(x) : ε(u(x)), div (σ(x)) = 0, ε = E, (1)
where . denot es t he spat ial average over V . T he boundary condit ions applied on ∂V should reproduce as closely as possible t he in sit u st at e of t he r.v.e.

T his st at e is very seldomly known and periodicit y condit ions are assumed t o close t he problem : t he local st rain field is decomposed int o it s average and a periodic fl uct uat ion ε(u(x)) = E + ε(u * (x)) where u * is periodic (not at ion u * #). T he t ensions σ.n t ake opposit e values on oppposit e sides of t he r.v.e.

(not at ion σ.n -#). T he reader is referred t o Suquet [START_REF] Suquet | Element s of Homogenizat ion for Inelast ic Solid Mechanics[END_REF] for more det ails about periodicit y condit ions.

T he local problem (1) can be re-writ t en by int roducing a homogeneous reference mat erial wit h elast ic st iffness c 0 :

σ(x) = c 0 : ε(u * (x)) + τ (x), div (σ(x)) = 0, u * #, σ.n -#, (2)
where t he polarizat ion field τ (x) is given by :

τ (x) = δc(x) : (ε(u * (x)) + E) + c 0 : E, δc(x) = c(x) -c 0 . (3)
T he solut ion of problem (2) can be expressed by means of t he periodic Green's operat or associat ed wit h t he elast icit y t ensor c 0 and reads, in real space and Fourier space respect ively :

ε(u * (x)) = -Γ 0 * τ (x), ε * (ξ) = - Γ0 (ξ) : τ (ξ) ∀ξ = 0, ε * (0) = 0.
Aft er subst it ut ing back t he expression (3) of τ in t his relat ion, t he init ial local problem (1) reduces t o t he periodic Lippm ann-Schwinger int egral equat ion for ε(u) which reads, in real space:

ε(u(x)) = -Γ 0 * (δc(x) : ε(u(x))) + E. (4)
T he operat or Γ 0 is explicit ely known in Fourier space for arbit rary anisot ropy of t he reference medium:

Γ0 ijkh = N 0 ik ξ j ξ h | (ijkh) , K 0 ik = c 0 ijkh ξ j ξ h , , N 0 = K 0 -1 , (5)
where t he symbol (ijkh) denot es symmet rizat ion wit h respect t o t he four indices (ijkh) result ing in minor and major symmet ries for Γ0 . Specific expressions for Γ 0 can be found in Mura [START_REF] Mura | Microm echanics of Defect s in Solids[END_REF] for different classes of anisot ropy of t he reference medium.

T he int egral equat ion ( 4) is solved by a fixed-point met hod :

ε(u i+1 ) = -Γ 0 * δc : ε(u i ) + E. (6)
T his algorit hm can be furt her simplified by not ing t hat Γ 0 * (c 0 : ε) = ε(u * ) and t he it erat ive met hod reads

ε(u i+1 ) = ε(u i ) -Γ 0 * σ i , σ i (x) = c(x) : ε(u i (x)). (7)
T he i + 1-t h it erat e of t he numerical algorit hm t ypically reads :

ε i and σ i being known a) σi = FT(σ i ), b) Convergence t est , c) εi+1 (ξ) = εi (ξ) - Γ0 (ξ) : σi (ξ) ∀ξ = 0 and εi+1 (0) = E, d) ε i+1 = FT -1 (ε i+1 ) e) σ i+1 (x) = c(x) : ε i+1 (x), ∀ x ∈ V.                (8)
FT st ands for t he Fourier t ransform. Convergence is reached when σ i+1 is in equilibrium.

In pract ice t he microst ruct ure is given in t he form of an image, consist ing of pixels (or voxels in dimension t hree) of a given size. T he spat ial resolut ion of t he image is t he number of pixels along each coordinat e axis. T he image sampling in real space generat es a corresponding sampling in Fourier space. T he algorit hm ( 8) is implement ed in discret ized form , t he FT being replaced by t he Fast Fourier Transform (FFT ) at t ached t o t he above sampling in real and Fourier space ( [START_REF] Moulinec | A numerical met hod for comput ing t he overall response of nonlinear composit es wit h complex microst ruct ure[END_REF]).

T he FFT met hod can be used in most problems where t he FEM is used, provided t hat periodicit y condit ions can be adopt ed for t he problem at hand. Ext ensions of t he met hod t o plast icit y wit h or wit hout hardening, t o viscoplast icit y, t o phase t ransformat ion were given in [START_REF] Moulinec | A numerical met hod for comput ing t he overall response of nonlinear composit es wit h complex microst ruct ure[END_REF] [7] [START_REF] Dreyer | A st udy of t he coarsening in P b/ Sn solders[END_REF] in t he cont ext of infinit esimal st rains. T he next sect ions are concerned wit h t he ext ension of t he met hod t o large st rains.

3.

Hyperelastic materials

Consider now t he case where t he deformat ions applied t o t he composit e mat erial are large. T he shape of t he r.v.e. V evolves wit h t ime and classically t wo configurat ions play a part icular role, t he init ial configurat ion V 0 and t he current configurat ion V (t). When t he individual const it uent s are hyperelast ic a Lagrangian formulat ion can be adopt ed and t he local problem can be posed on t he init ial configurat ion V 0 only. A part icle which was init ially at point X moves t o a locat ion x = x(X, t) = X + u(X, t) at t ime t. T he deformat ion gradient and t he first P iola-Kirchhoff st ress t ensor π are relat ed by:

f = I + ∇ X u(X, t), π(X) = ∂w ∂f (X, f ), (9)
where w(X, .) is t he st rain energy at point X. T he average gradient F of t he t ransformat ion is applied increment ally along a prescribed pat h in t he loading space. T his pat h is paramet rized by a scalar variable t. T he local problem t o be solved reads :

π(X) = ∂w ∂f (X, f ), div X π = 0 in V 0 , f (X) = F + ∇ X u * , u * #, π.N -# on ∂V 0      (10)
where N denot es t he out er normal unit vect or t o ∂V 0 . T he local problem (10) is similar t o (2). It is posed on a fixed configurat ion which allows us t o use t he FFT met hod wit h fixed grids in real space and Fourier space. T he const it ut ive relat ions are nonlinear and t he whole gradient f (and not only it s symmet ric part ) are t aken int o account int o [START_REF] Lebensohn | N-sit e modelling of a 3d viscoplast ic polycryst al using Fast Fourier T ransforms[END_REF]. T he problem [START_REF] Lebensohn | N-sit e modelling of a 3d viscoplast ic polycryst al using Fast Fourier T ransforms[END_REF] is solved st ep-by-st ep in t ime, wit h a Newt on-Raphson algorit hm at each t ime st ep.

Newton-Raphson algorithm

F t+∆t is imposed at t ime t + ∆t. T he principal unknown is t he gradient ∇u * t+∆t of t he periodic fl uct uat ion of t he displacement . T his gradient is det ermined by imposing t hat π t+∆t is in equilibrium.

Let π i , f i and (u * ) i denot e t he it erat es approximat ing π t+∆t , f t+∆t and u * t+∆t .

It erat e i: ∇(u * ) i-1 being known :

(1) Comput e f i-1 and π i-

1 f i-1 (X) = F t+∆t + ∇(u * ) i-1 (X), π i-1 (X) = ∂w ∂f (X, f i-1 (X)). ( 11
)
Check if π i-1 is in equilibrium. If not , t hen

(2) Solve t he linear t angent problem for t he periodic field ∇(δu * )

∂ 2 w ∂f 2 (f i-1 ) : ∇(δu * ) : ∇v = -π i-1 : ∇v ∀v #. (12)
(3) Updat e ∇(u

* ) i : ∇(u * ) i (X) = ∇(u * ) i-1 (X) + ∇(δu * )(X).
T he crit erion which serves t o check equilibrium reads :

max ξ |ξ. π i (ξ)| ≤ ε| π i (0)|, where t ypicallyε = 10 -4 .
T he problem ( 12) is solved it erat ively using t he FFT met hod described in sect ion 2. More specifically an int ernal loop is performed t o find ∇(δu * )(X)

∇(δu * ) k = ∇(δu * ) k-1 - Γ0 * L i-1 : ∇(δu * ) k-1 + π i-1 , (13)
where

L i-1 (X) = ∂ 2 w ∂f 2 (X, f i-1 (X)).
Γ0 is defined in Fourier space as

Γ0 ijkh = N 0 ik ξ j ξ h | (ik)(jh) ,
where t he symbol (ik)(jh) denot es symmet rizat ion wit h respect t o t he indices i, k and j, h only.

Example: fiber-reinforced elastomer

T he above algorit hm has been applied t o model t he deformat ion of an elast omeric mat rix reinforced by st iff fibers. In it s init ial configurat ion t he unit cell was a square cont aining 64 circular ident ical impenet rable fibers, wit h 25% volume fract ion, arranged randomly in t he unit cell. T he mat rix and t he fibers were compressible Mooney-Rivlin mat erials wit h st rain energy :

w(f ) = µ 2 (i 1 -3) -µ ln j + κ 2 (j -1) 2 , i 1 = tr( T f .f ), j = det f . ( 14
)
T he fibers were 10 t imes st iffer t han t he mat rix. 10 different configurat ions were t est ed. T he macroscopic deformat ion was a biaxial isochoric deformat ion F = exp (E 1 ) e 1 ⊗ e 1 + exp (-E 1 ) e 2 ⊗ e 2 + e 3 ⊗ e 3 , where t he logarit hmic st rain in t he first direct ion E 1 was negat ive (cont ract ion in direct ion 1, ext ension in direct ion 2). T he deformed st at es of a t ypical configurat ion are shown in figure 1. T he st ress-st rain response of all 10 different configurat ions are shown in Figure E1 = Macroscopic st ress-st rain curves. First P iola-Kirchhoff st ress as a funct ion of t he logarit hmic st rain in direct ion 1. 10 different configurat ions. Comparison bet ween t he FFT simulat ions and t he second order est imat e (SOE) wit h t he lower Hashin-Sht rikman est imat e.

2. T hey are seen t o be very close t o one anot her. Also shown in Figure 2 is t he predict ion of t he second-order est imat e (SOE) of P ont e Cast a ñeda [START_REF] Castañeda | A second-order homogenizat ion met hod in finit e elast icit y and applicat ions t o black-filled elast omers[END_REF] implement ed wit h t he Hashin-Sht rikman est imat e (lower bound) for circular fibers in a mat rix. As can be seen from t his figure, t he agreement bet ween t he t heory and t he numerical calculat ions is excellent .

4.

Eulerian formulation

Local problem

We now t urn t o t he case where t he const it ut ive relat ions of t he const it uent s can be expressed by means of quant it ies defined on t he current configurat ion only, t he Cauchy st ress σ and t he Eulerian st rain-rat e d for inst ance. T his is t he case for Newt onian fl uids for which t he const it ut ive relat ions read

σ(x) = L(x) : d(x), (15)
and more generally for viscous mat erials when elast icit y effect s are neglect ed (t he case of power-law viscous mat erials is considered in Lebensohn [START_REF] Lebensohn | N-sit e modelling of a 3d viscoplast ic polycryst al using Fast Fourier T ransforms[END_REF]). For simplicit y we rest rict our at t ent ion t o Newt onian fl uids as described by (15). T he geomet ry of t he r.v.e. evolves wit h t ime in t wo ways. First t he unit vect ors Y i defining t he unit cell are convect ed by t he macroscopic velocit y gradient :

Ẏ i (t) = grad v(t) .Y i (t). ( 16
)
Second, t he posit ion of t he phases wit hin t he r.v.e. changes wit h t ime, t he mat erial part icles being convect ed by t he microscopic velocit y field:

ẋ(t) = v(x, t) (17)
T he local problem t o be solved in t he Eulerian configurat ion consist s of t he const it ut ive relat ions, t oget her wit h t he equilibrium and compat ibilit y equat ions. At each t ime t t he local problem on t he current configurat ion can be solved using t he FFT met hod described in sect ion 2. It remains t o see how t he unit cell changes wit h t ime.

A Multi-Particle Method

The two grids. T he proposed met hod falls in t he cat egory of mult i-part icle met hods and follows rat her closely t he P art icle-In-Cell (P IC) met hod proposed by Sulsky, Chen and Schreyer [START_REF] Sulsky | A part icle met hod for hist ory-dependent mat erials. Com p. Met hods[END_REF][START_REF] Sulsky | Applicat ion of a part icle-in-cell met hods t o solid mechanics[END_REF], wit h a difference st emming from t he updat ing of t he comput at ional grid (t he met hod proposed by Lebensohn [10] is also similar in spirit ). T he underlying idea of t he P IC met hod is t o consider t wo separat e grids:

-T he com put at ional grid is used for applying t he FFT met hod described in sect ion 2. It is a regular grid but it does not have t o be rect angular (Fourier t ransforms can be defined on noncubic lat t ices).

-T he m at erial grid is at t ached t o t he mat erial part icles. It does not have t o be st ruct ured or regular and can be seen as a collect ion of part icles, rat her t han as a st ruct ured mesh. It is only used t o apply t he const it ut ive relat ions (which are mat erial relat ions). In t he FFT met hod it is used in real space only. In most met hods t he mat erial grid is finer t han t he comput at ional one (by a fact or of 4 t o 9).

Each grid carries it s own set of unknowns. In ot her words, -t he comput at ional grid carries unknowns labelled wit h lowerscript c (as in " comput at ions" ) :

x (i,j) c = (i -1)Y 1 + (j -1)
Y 2 denot e t he nodes of t he comput at ional grid, σ c , v c denot e t he st ress and velocit y fields at t hese nodes.

-t he mat erial grid carries unknowns labelled wit h lowerscript p (as in " part icles" ): x p denot e t he locat ion of t he part icles, σ p , v p denot e t he st ress and velocit y field at t he part icles.

Updating the grids

In t he init ial P IC met hod ( [START_REF] Sulsky | A part icle met hod for hist ory-dependent mat erials. Com p. Met hods[END_REF][START_REF] Sulsky | Applicat ion of a part icle-in-cell met hods t o solid mechanics[END_REF]), t he comput at ional grid is fixed. However in t he problem under considerat ion here t he comput at ional grid is associat ed wit h t he unit cell which generat es t he whole microst ruct ure by periodicit y and cannot remain fixed. T he unit vect ors Y i defining t he periodic lat t ice are mat erial vect ors and follow t he m acroscopic velocit y field according t o (16). T he comput at ional grid is updat ed accordingly. T he mat erial grid is updat ed using t he m icrosocopic velocit y field according t o (17).

T he grids are updat ed using an explicit scheme. At t ime t n , t he velocit y fields v c (t n ) is comput ed on t he comput at ional grid and t ransferred (see below) ont o t he mat erial grid t o get a field v p (t n ). T hen :

Y i (t n+1 ) = Y i (t n ) + ∆t grad v c (t n ) .Y i (t), x c (t n+1 ) = (i -1)Y 1 (t n+1 ) + (j -1)Y 2 (t n+1 ), x p (t n+1 ) = x p (t n ) + ∆t v p (x p (t n )).    ( 18
)
Transfer operators. Transfer operat ors are needed in order t o t ransfer informat ion from t he comput at ional grid t o t he mat erial grid and vice-versa. Devising consist ent t ransfer operat ors is cert ainly a crucial part of t he algorit hm.

Consider a comput at ional grid wit h an init ially rect angular pat t ern. In t he init ial configurat ion t he grid nodes have coordinat es ((i -1)∆x 0 , (j -1)∆y 0 ). In t he current configurat ion t he grid nodes have coordinat es (x c i,j , y c i,j ). Let us int roduce a family of shape funct ions N (i,j) , i = 1, ..., I, j = 1, ..., J (as is classical in t he FEM) on t he comput at ional grid, where I and J denot e t he number of sampling point s on t he horizont al and vert ical axis (for t he init ial grid) respect ively. Given t he regular geomet ry of t he comput at ional grid a nat ural choice for t he N (i,j) are t he int erpolat ion funct ions of t he 4 point quadrilat eral element in t he FEM (t hese funct ions were used in t he example present ed in sect ion 4.4). T hen, given t he values f (i,j) c = f (x c (i, j) of any funct ion f at t he comput at ional nodes, t he int erpolat ed field at any ot her point x read as

f (x) = I,J (i,j)=(1,1) N (i,j) (x)f (i,j) c . (19)
T hese relat ions apply in part icular at t he part icle posit ions x p . T he inverse t ransfer from t he mat erial grid t o t he comput at ional grid is performed using t he same int erpolat ion funct ions. Each part icle is assigned a mass m p and t he mass of a comput at ional node is defined as

m(x (i,j) c ) = P p=1 m p N (i,j) (x p ) (20)
T hen t he values f c of f at any comput at ional node x c is given by :

f c (x (i,j) c ) = 1 m(x (i,j) c ) P p=1 m p N (i,j) (x p )f p (x p ) (21)

A test example

To check t he accuracy of t he updat ing scheme, t he problem of a rigid part icle rot at ing in a shear fl ow of a newt onian fl uid wit h uniform and const ant viscosit y has been invest igat ed. T he fl uid is subject ed t o a macroscopic fl ow wit h st rainrat e D and rot at ion-rat e Ω. Inert ia effect s are neglect ed. T he solut ion (rot at ion of t he part icle) for a spheroidal part icle wit h infinit esimal volume fract ion c (p) in an incompressible newt onian fl uid goes back t o J effery [START_REF] Jeffery | T he mot ion of ellipsoidal part icles immersed in a viscous fl uid[END_REF]. A generalizat ion of J effery's result , also applying t o compressible fl uids and viscous part icles, can also be derived by means of Eshelby's result for ellipsoidal inclusions. T his derivat ion will not be given here. For a rigid part icle t he evolut ion equat ion for a unit normal vect or u at t ached t o t he part icle is

du dt = ω I .u, ω I = Ω -Π : S -1 : D, (22)
where S and Π are t he t wo Eshelby t ensors of t he part icle in t he mat rix giving respect ively t he deformat ion-rat e in t he inclusion and t o t he rot at ion-rat e of t he inclusion.

T his general t hree-dimensional solut ion can be specialized t o dimension 2. T he result is

du dt = Ω.u + B (D.u -(u.D.u)u) , (23)
wit h B = (r 2 -1)/ r 2 + 1 + 2r 1-2ν 2(1-ν) where r is t he aspect -rat io of t he part icle and ν is t he P oisson rat io of t he linearly viscous mat rix. In a shear fl ow,

x 1 (t) = X 1 + γ(t)X 2 , x 2 (t) = X 2 , (24)
(23) gives t he t he rot at ion of t he unit normal vect or u = sin ϕe 1 + cos ϕe 2 in t he form :

dϕ dt = γ 2 1 + B(cos 2 ϕ -sin 2 ϕ) , (25)
which can be int egrat ed int o

ϕ(t) = ϕ(0) + Arct g R t g γ(t) R + 1/R , wit h R = 1 + B 1 -B . (26)
A numerical simulat ion of t he rot at ion of a part icle in a shear fl ow has been conduct ed using t he met hod described in sect ion 4.2 wit h t he following dat a r = 3, ν = 0.499, c (p) = 0.3%.

T he comparison bet ween t he analyt ical and t he numerical result s is good as shown in Figure 4 b. Also shown in t his figure is t he lat t ice rot at ion. As can be seen t he part icle rot at es fast er t han t he lat t ice.

Concluding remarks

Two ext ensions t o finit e st rains of t he numerical met hod based on Fast Fourier Transforms ( [START_REF] Moulinec | A fast numerical met hod for comput ing t he linear and nonlinear propert ies of composit es[END_REF][START_REF] Moulinec | A numerical met hod for comput ing t he overall response of nonlinear composit es wit h complex microst ruct ure[END_REF]) have been proposed and t est examples have been analyzed. A few numerical issues deserve at t ent ion in fut ure work:

In t he Lagrangian approach present ed in sect ion 3, t he rat e of convergence of t he it erat ive FFT algorit hm applied t o t he t angent problem (13) can be rat her poor. T his is probably due t o t he fact t hat t he t angent operat or is only st rongly ellipt ic (and not very st rongly ellipt ic) and highly cont rast ed. An accelerat ed scheme would be very helpful.

In t he Eulerian approach present ed in sect ion 4, consist ent t ransfer operat ors bet ween t he comput at ional grid and t he mat erial grid are t o be found.

Figure 1 .

 1 Figure 1. Hyperelast ic mat rix reinforced by 64 circular fibers. (a): init ial configurat ion. (b)-(f): successive deformed st at es for different values of t he macroscopic logarit hmic st rain E1 = ln F11.

Figure 2 .

 2 Figure 2.Macroscopic st ress-st rain curves. First P iola-Kirchhoff st ress as a funct ion of t he logarit hmic st rain in direct ion 1. 10 different configurat ions. Comparison bet ween t he FFT simulat ions and t he second order est imat e (SOE) wit h t he lower Hashin-Sht rikman est imat e.

Figure 3 .

 3 Figure 3.Top row: init ial grids. Left : t he com put at ional grid. Cent er: t he m at erial grid. R ight : bot h grids superim posed. Bot t om row: convect ed grids.

Figure 4 .

 4 Figure 4. Part icle rot at ing in a shear fl ow.

Acknowledgments

T he aut hors great ly benefit ed from discussions wit h R. Lebensohn, D. Sulsky and N. Triant afyllidis during t he complet ion of t his research. P art of it was carried out while P ierre Suquet was a Visit ing P rofessor at Calt ech. Wonderful int eract ions wit h K. Bhat t acharya and M. Ort iz, financial support from Calt ech and from t he Délégat ion G én érale à l'A rm em ent are grat efully acknowledged.


St Weijun Wang
email: weijun.wang@irccyn.ec-nantes.fr

Fouad Caro

Bennis

Oscar Roberto
email: oscarrobertosalinas@hotmail.com

Salinas Mejia

A Simplified Morphing Blade for Horizontal Axis Wind Turbines

The aim of designing wind turbine blades is to improve the power capture ability. Since rotor control technology is currently limited to controlling rotational speed and blade pitch, an increasing concern has been given to morphing blades. In this paper, a simplified morphing blade is introduced, which has a linear twist distribution along the span and a shape that can be controlled by adjusting the twist of the blade's root and tip. To evaluate the performance of wind turbine blades, a numerical code based on the blade element momentum theory is developed and validated. The blade of the NREL Phase VI wind turbine is taken as a reference blade and has a fixed pitch. The optimization problems associated with the control of the morphing blade and a blade with pitch control are formulated. The optimal results show that the morphing blade gives better results than the blade with pitch control in terms of produced power. Under the assumption that at a given site, the annual average wind speed is known and the wind speed follows a Rayleigh distribution, the annual energy production of wind turbines was evaluated for three types of blade, namely, morphing blade, blade with pitch control and fixed pitch blade. For an annual average wind speed varying between 5 m/s and 15 m/s, it turns out that the annual energy production of the wind turbine containing morphing blades is 24.5 % to 69.7 % higher than the annual energy production of the wind turbine containing pitch fixed blades. Likewise, the annual energy production of the wind turbine containing blades with pitch control is 22.7 % to 66.9 % higher than the annual energy production of the wind turbine containing pitch fixed blades. Wind energy is growing more and more popular worldwide. The search for ways to make wind turbines (WT) more efficient and competitive becomes paramount. The main control methods used to optimize or limit the power extracted from the wind turbine are usually based on the control of the rotor's rotational speed and the pitch of the blades. Recently, an increasing concern has been given to Morphing Blades (MB) [START_REF] Lobitz | Enhanced performance of HAWTs using adaptive blades[END_REF][START_REF] Beyene | A Morphing Blade for Wave and Wind Energy Conversion[END_REF][START_REF] Daynes | A Morphing Wind Turbine Blade Control Surface[END_REF].

It is well known that MB can improve the power efficiency by changing their shape according to variations in wind speed. Moreover, they have the potential to significantly relieve unwanted stresses in the blades to prolong their life while the wind is very harsh on them. Barlas et al. [4] presented status of active aero-elastic rotor control research for wind turbines in terms of using advanced control concepts to reduce loads on the rotor. They analyzed the smart control concepts, including twist control, camber control and moveable control surfaces (trailing edge flaps or servo tabs actuated by smart materials). In the scope of this study, the term "morphing blade" means that the angle of attack of each section of the blade is controlled rather than the aerodynamic characteristic curve of the section.

MBs can be either "passive" or "active" depending on the type of twist control. The former are essentially monolithic and rely on the flexibility of their structure or, in other words, on elastic deformation. In contrast, "active adaptive blades" are made up of a number of independent span-wise sections that can be oriented to achieve any desired twist distribution. Although the adaptation of the blade shape has the ability to optimize the efficiency, due to the complex structure and high cost, its application has been limited. In this paper, a simplified MB is introduced. This morphing blade has a linear twist distribution along the span.

Since the combined Blade Element Momentum Theory (BEMT) is a fairly accurate analytical tool and has low computational cost, it is widely used in the wind energy industry to estimate the theoretical output power from a rotor with defined blade dimensions [START_REF] Duran | Computer aided design of horizontal axis 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] Kulunk | HAWT Rotor Design and Performance Analysis[END_REF][START_REF] Tenguria | Investigation of blade performance of horizontal axis wind turbine based on blade element momentum theory ( BEMT ) using NACA airfoils[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]. A numerical code is developed in this paper based on the BEMT and validated by comparison with the experimental results of the NREL phase VI test turbine [START_REF] Hand | Unsteady Aerodynamics Experiment Phase V: Test Configuration and Available Data Campaigns[END_REF][START_REF] Lanzafame | Horizontal axis wind turbine working at maximum power coefficient continuously[END_REF], which is a stall-controlled wind turbine.

This turbine is regulated using passive stall methods at high wind speeds to limit the output power. Thus, there is great potential for improvement of its power capture capability. Here, we use the blade of the NREL Phase VI WT as a reference blade, and formulate the optimization problem associated with the proposed simplified MB. To make a good comparison, the optimization problem associated with the blade with pitch control (BPC) is also formulated. The optimal results show that the MB gives better results than the BPC in terms of produced power. Under the assumption that in a given site, the annual average wind speed is known and the wind speed follows a Rayleigh distribution, we can evaluate the annual energy produced by these three types of blade.

The paper is organized as follows. The BEMT is presented in Sec. 2. The calculation model used to evaluate the WT produced power is described in Sec. 3. The proposed simplified MB is introduced in Sec. 4. Some optimization problems are formulated in Sec. 5 for the control of the wind turbines. Optimal results are also given. Finally, some concluding remarks and future work are provided in Sec. 6.

Blade Element Momentum Theory 2.1 Momentum Theory

Fig. 1: Schematic of momentum theory for wind turbines Figure 1 gives a schematic of the momentum theory for wind turbines. Based on some assumptions, a simple model, known as actuator disc model, can be used to determine the power from an ideal turbine rotor and the thrust of the wind on the ideal rotor [START_REF] Duran | Computer aided design of horizontal axis wind turbine blades[END_REF][START_REF] Kulunk | HAWT Rotor Design and Performance Analysis[END_REF]. From axial momentum and angular momentum, the element of thrust dF a and the element of torque dT t can be obtained as:

dF a = 4πρv 2 0 a(1 -a)rdr (1)
dT t = 4πρv 0 ωa ′ (1 -a)r 3 dr (2)

Blade Element Theory

In order to apply blade element analysis, it is assumed that the blade is divided into N sections. The analysis is based on some assumptions, such as that there is no aerodynamic interaction between different blade elements and the forces on the blade elements are solely determined by the lift and drag coefficients [START_REF] Duran | Computer aided design of horizontal axis wind turbine blades[END_REF][START_REF] Kulunk | HAWT Rotor Design and Performance Analysis[END_REF][START_REF] Manwell | Wind Energy Explained Theory,Design and Application[END_REF].

Figure 2 shows the velocities and forces on a blade element of a wind turbine blade. As a result, the following equations are obtained:

φ = arctan v 0 (1 -a) ωr(1 + a ′ ) (3)
v r = v 0 (1 -a) sin φ (4) α = φ -γ (5)
Lift and drag forces on an annular blade element are given by:

dF L = ρc 2 v 2 r C L dr (6)
dF D = ρc 2 v 2 r C D dr (7)
These lift coefficient C L and drag coefficient C D depend on the angle of attack α and the blade profile. The element of thrust dF a and the element of torque dT t are expressed as:

dF a = b ρc 2 v 2 r (C L cos φ + C D sin φ)dr ( 8
)
dT t = rdF t = b ρc 2 v 2 r (C L sin φ -C D cos φ)rdr (9)
4πρv 2 0 a(1 -a)rdr = b ρc 2 v 2 r (C L cos φ + C D sin φ)dr(10) 4πρv 0 ωa ′ (1 -a)r 3 dr = b ρc 2 v 2 r (C L sin φ -C D cos φ)dr(11)
After some algebraic manipulations and by adding the correction of the Prandtl tip loss factor, the following relationships are obtained [START_REF] Lanzafame | Fluid dynamics wind turbine design: Critical analysis, optimization and application of BEM theory[END_REF]:

a = 1 4F sin 2 φ σ(C L cosφ + C D sin φ) + 1 ( 12
)
a ′ = 1 4F sin φ cos φ σ(C L sin φ -C D cos φ) + 1 ( 13
)
where F is the Prandtl tip loss factor defined as:

F = 1 π arccos[exp( b(r -r t ) 2r sin φ )] (14)
and σ is the rotor solidity, defined as:

σ = cb 2πr ( 15
)
It is to be noted that Eqn. ( 12) is valid for the axial induction factor value to be between 0 and 0.4. For axial induction factor greater than 0.4, there are several methods to obtain it [START_REF] Lanzafame | Fluid dynamics wind turbine design: Critical analysis, optimization and application of BEM theory[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][START_REF] Buhl | A New Empirical Relationship between Thrust Coefficient and Induction Factor for the Turbulent Windmill State[END_REF][START_REF] Lanzafame | Design and performance of a double-pitch wind turbine with non-twisted blades[END_REF]. In this paper, we use the method proposed by Buhl [START_REF] Buhl | A New Empirical Relationship between Thrust Coefficient and Induction Factor for the Turbulent Windmill State[END_REF] and implemented by R. Lanzafame [START_REF] Lanzafame | Fluid dynamics wind turbine design: Critical analysis, optimization and application of BEM theory[END_REF], namely

a = 18F -20 -3 C T (50 -36F) + 12F(3F -4) 36F -50 ( 16
)
where C T is the thrust coefficient of the wind turbine rotor.

For each blade element, it can be calculated as [START_REF] Duran | Computer aided design of horizontal axis wind turbine blades[END_REF]:

C T = dF a 1/2ρv 2 0 2πrdr = σ(1 -a) 2 (C L cos φ + C D sin φ) sin 2 φ ( 17
)
3 Calculation Model 3.1 Model Formulation One of the most difficult issues for BEMT is to determine the induction factors (a and a ′ ) and the correct lift and drag coefficients (C L and C D ).

Here, we determine the induction factors by using an iteration method, which is currently the governing method [START_REF] Duran | Computer aided design of horizontal axis wind turbine blades[END_REF][START_REF] Kulunk | HAWT Rotor Design and Performance Analysis[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]. Figure 3 shows the flowchart to determine the induction factors a and a ′ . The parameters associated with each section are given, namely, the chord length c, the twist angle γ, the free stream wind speed v 0 , the air density ρ, the radial coordinate at rotor plane r, the root radius r t of the wind turbine, the rotor rotational speed ω, the number of blades b, they airfoil type and its aerodynamic parameters.

Fig. 3: Calculation flowchart for induction factors

In order to find a piecewise polynomial relation between α and C D (α and C L , respectively) from the experimental data, R. Lanzafame implemented a fifth-order logarithmic polynomial for the angle of attack from -6 • to 20 • . For the angle of attack from 20 • to 45 • , the mathematical functions of the following equations were implemented [START_REF] Lanzafame | Fluid dynamics wind turbine design: Critical analysis, optimization and application of BEM theory[END_REF]:

C L = 2C L max sin α cos α (18)
C D = 2C D max sin 2 α (19)
The C L max and C D max are shown as follows:

C L max = C L | α=45 • and C D max = C D | α=90 • (20)
In this work, a cubic spline interpolation was implemented for the angle of attack from -20 • to 20 • , and the mathematical functions in Eqn. [START_REF] Vieira Da Rosa | Fundamentals of Renewable Energy Processes[END_REF] and Eqn. (19) were implemented for the angle of attack from 20 • to 90 • .

The experimental data for a two-dimensional S809 airfoil section were obtained obtained at the OSU (Ohio State University) wind tunnel with a Reynolds number from 990,000 to 1,040,000 [START_REF] Hand | Unsteady Aerodynamics Experiment Phase V: Test Configuration and Available Data Campaigns[END_REF][START_REF] Ramsay | Effects of Grit Roughness and Pitch Oscillations on the S809 Airfoil: Airfoil Performance Report[END_REF][START_REF] Ramsay | Wind Tunnel Testing of Three S809 Aileron Configurations for use on Horizontal Axis Wind Turbines[END_REF]. Figure 4 shows C D as a function of α and C L as a function of α. This figure also provides a comparison of the functions from the fifth-order logarithmic polynomial and that from the cubic spline interpolation.

F a = b ∑ j=1 r t r r dF a ( 21
)
T t = b ∑ j=1 r t r r dT t (22)
Then the produced power P and the power coefficient C p of a wind turbine rotor are expressed as follows:

P = T t ω (23)
C p = P 1 2 ρπr 2 t v 3 0 (24)

Validation of the Model

To verify the validity of the formulated calculation model, we take the NERL Phase VI WT as an example. This wind turbine has two twisted blades, a variable chord along the blade, and a rotor of diameter equal to 10 m. The blade section is of type S809. The rotational speed of the rotor is constant and equal to 72 rpm with a rated produced power equal to 19.8 kW [START_REF] Hand | Unsteady Aerodynamics Experiment Phase V: Test Configuration and Available Data Campaigns[END_REF][START_REF] Lanzafame | Horizontal axis wind turbine working at maximum power coefficient continuously[END_REF].

In accordance with experiments, the angle between the rotor plane and the tip chord is constant and equal to 3 deg. Therefore, we end up with a Fixed Pitch Blade (FPB) and the fixed pitch angle is equal to 4.775 deg. Figure 5 shows the distribution of the actual twist and chord distribution of the FPB.

In the calculation model, the root radius r r and tip radius r t of the blade are equal to 1.27 m and 5 m, respectively. The blade is divided into 18 cross-sections. The air density ρ is equal to 1.25 kg/m 3 and the reference wind speed varies 5 m/s and 25 m/s. In order to make a better comparison, in Fig. 6 we also provide the obtained simulation results of the calculation model representing the lift and drag coefficients based on a fifth-order logarithmic polynomial (the angle of attack from -20 • to 20 • ). Figure 6 shows a good agreement between simulated and experimental results. Usually, the well controlled WTs operate in two primary regimes [START_REF] Merabet | Torque and Pitch Angle Control for Variable Speed Wind Turbines in All Operating Regimes[END_REF], namely, the partial load region and the full load region, as shown in Fig. 7. The cut-in wind speed v ci is the minimum wind speed at which the WT generates usable power and the cut-out wind speed v co is the one at which the braking system is activated to slow down or stop the turbine in order to avoid any damage.

Meanwhile, to protect the components of the WT, wind power must be shed to limit output power when wind speed is over the rated speed. The partial load region where wind speed varies between cut-in speed v ci and rated speed v rated is named "Region 1". The full load region where wind speed varies between rated speed v rated and cut-out speed v co is named "Region 2".

The Fixed Pitch Blade

It is well known that due to geometric reasons, the effective wind velocity vector varies both, in magnitude and direction, along a wind turbine blade. An optimum performance -assuming constant airfoil section throughout the blade -would require the exact same relative orientation between the local wind velocity vector and the corresponding cross-section of the blade. As an attempt to reach s-uch optimum performance, modern blades are produced with a physical twist, which may correspond to the theoretically ideal value.

A Fixed Pitch Blade (FPB) operates at a fixed pitch angle. For a FPB, the actual twist angle γ i of the ith blade ele-Fig. 7: Produced power as a function of wind speed ment is the sum of its pre-twist angle γ i0 and the fixed pitch angle β 0 , namely,

γ i = γ i0 + β 0 (25)

The Blade with Pitch Control

A Blade with Pitch Control (BPC) implies that the pitch of the WT's rotor blades can be adjusted by a pitch control system. An optimal design is needed to pitch the rotor blades in order to maximize the produced power at any wind speed. In Region 1, pitch control can improve the power coefficient. In Region 2, pitch control can maintain a constant output power, P rated . Fig. 8: Schematic of a blade with pitch control Figure 8 shows the schematic of a BPC. Since the pitch control system can not change the shape of the blades, any adjustment will result in all blades sections experiencing the same change in actual twist angle. For a twisted blade, the actual twist angle of the ith element is the sum of the pretwist angle and the pitch control angle, β:

γ i = γ i0 + β (26)

The Simplified Morphing Blades

A morphing blade has the ability to approach the optimal profile and especially the optimal angle of attack for each element. This aim can be achieved in an active or passive manner.

Figure 9 shows a MB with active control. This blade has a constant airfoil section throughout the blade and the twist angle of each element can be controlled by some mechanisms, to correspond to the theoretically ideal value. It is apparent that the mechanisms used to control the elements individually may be complex and expensive.

γ i = (γ r -γ t ) r i -r t r r -r t + γ t ( 27
)
It is noteworthy that the conventional blade has a fixed shape and a hollow profile usually formed by two shell structures and some webs, as shown in Fig. 11. For MB, since the shape of the blade is changeable, the blade must be divided into several elements and there must be a flexible skin along the span.

Optimization Problems and Result Analysis 5.1 Optimization Problem for the Simplified Morphing Blade

The Phase VI WT is taken as a reference [START_REF] Hand | Unsteady Aerodynamics Experiment Phase V: Test Configuration and Available Data Campaigns[END_REF][START_REF] Lanzafame | Horizontal axis wind turbine working at maximum power coefficient continuously[END_REF]. The simplified morphing blade is divided into 18 cross-sections. The objective is to maximize the output power P at a steady wind speed. For the simplified morphing blade, the two twist angles at the first element and at the last element, γ r and γ t are decision variables. Then the actual twist angle of the ith blade element can be calculated as Eqn. p is the set of design parameters given in Table 1. The only constraint is that the output power should be less than or equal to the rated power : P ≤ 19.8 kW .

As a consequence, for a steady reference wind speed v re , which varies from v ci to v co , the optimization problem is for- The performance function, the other design parameters and the constraint are the same as in the former optimization problem associated with the MB.

Then, for a steady reference wind speed, the optimization problem for the blade with pitch control is formulated as follows:

minimize -P(x, p) over x = β p = [γ γ γ b c r t r r ω ρ v re ] subject to P ≤ 19.8 kW -5deg ≤ β ≤ 25deg (29)

Result Analysis 5.3.1 Optimization Results

The MATLAB f mincon function was used to solve optimization problems (28) and (29). Several starting points were used to come up with results, which are as close as possible to global optima.

Figure 12 shows the optimal performance of the MB and BPC at different wind speeds. Moreover, performance of the fixed pitch turbine, the Phase VI WT with a fixed pitch, are given.

The results show that the BPC can improve the power coefficients C p when wind speed is higher than 9 m/s. The main reason is that the FPB is stalled to maintain the output power when wind speed is higher than 12 m/s . Moreover, the results show that the produced power is higher with the MB than with BPC, except for some wind speed ranges. Since the actual twist angles of the MB are linear and the actual twist angles of the BPC are non-linear, it is understandable that the performance of BPC is better than MB for some wind speeds.

Table 2 shows the maximum power produced by wind turbines containing FPB, MB and BPC for different wind speeds. Besides, the optimum twist angles γ r and γ t are given for the MB and the optimum pitch control angle β is given for the BPC for the different wind speeds. The fixed pitch angle β 0 for the FPB is equal to 4.775 deg.

Figure 13 depicts the actual twist angles as a function of the blade radius for the three types of blade under study

Comparing Annual Energy Production

The Annual Energy Production (AEP) can be used as an index for the comparison of wind turbines. The AEP depends on the power curve of the wind turbine at hand and the variability of the wind.

The wind speed can be modeled by using a Weibull Distribution. The probability curve is usually defined by two parameters: the average wind speed v and the shape factor k [START_REF] Vieira Da Rosa | Fundamentals of Renewable Energy Processes[END_REF]. A simple model is chosen to assess the wind speed frequency. Indeed, a Rayleigh distribution is considered by setting the shape factor k to 2, i.e., k = 2. Therefore, for a given average wind speed v, the probability p(v) for the wind speed to be equal to v is obtained by the following formula:

p(v) = π 2 v v 2 exp - π 4 v v 2 (30)
Figure 14 illustrates the Rayleigh distributions of the wind speed for three average wind speeds, i.e., v = 5 m/s, v = 10 m/s and v = 15 m/s. Table 3 and Fig. 15 show the annual energy production of the wind turbines as a function of the average wind speed and for the three types of blade. It appears that the AEP of the wind turbine containing morphing blades is 24.5 % to 69.7 % higher than the AEP of the wind turbine containing pitch fixed blades. Likewise, the AEP of the wind turbine contain- ing blades with pitch control is 22.7 % to 66.9 % higher than the AEP of the wind turbine containing pitch fixed blades.

Conclusions and Future Work

A simplified morphing blade was introduced in this paper in order to improve the power capture ability of wind turbines. This morphing blade has a linear twist distribution along the span and its shape can be controlled by adjusting the root twist angle and the tip twist angle of the blade. A simulation model to compute the performance of wind turbines based on blade element momentum theory has been formulated and verified. Then, the optimization problems to maximize the produced power, associated with the morphing blades and the blades with pitch control were both formulated. At different wind speeds, there must be different adjustable parameters for the morphing blades and the blades with pitch control, which are decision variables of the optimization problems at hand. The optimal results show that the morphing blades provides better results than the blades with pitch control in terms of produced power. Under the assumption that at a given site, the annual average wind speed is known and the wind speed follows a Rayleigh distribution, the annual energy production of wind turbines was evaluated for three types of blade, namely, morphing blade, blade with pitch control and fixed pitch blade. For an annual average wind speed varying between 5 m/s and 15 m/s, it turns out that the annual energy production of the wind turbine containing morphing blades is 24.5 % to 69.7 % higher than the annual energy production of the wind turbine containing pitch fixed blades. Likewise, the annual energy production of the wind turbine containing blades with pitch control is 22.7 % to 66.9 % higher than the annual energy production of the wind turbine containing pitch fixed blades. Finally, the detailed design of the inner structures of the proposed simplified morphing blade is part of the future work.

Nomenclature C p

 p Power coefficient of wind turbine rotor C T Thrust coefficient of wind turbine rotor v Wind speed v 0 Free stream wind speed v r Relative wind speed v re Reference wind speed v ci Cut-in wind speed v co Cut-out wind speed v rated Rated wind speed of wind turbine v Average wind speed p(v) Probability for the wind speed to be equal to v a Axial induction factor at rotor plane a ′ Angular induction factor b Number of blades of a rotor N Number of blade elements ρ Air density P Produced Power of wind turbine rotor P rated Rated Produced Power of wind turbine rotor x Vector of the decision variables p Vector of the design parameters r Radial coordinate at rotor plane r t Tip Radius of the blade r r Root radius of the blade r i Blade radius for the ith element F D Drag force on an annular blade element F L Lift force on an annular blade element F a The axial force on the blade element F t The edgewise forces on the blade element T t Rotor torque C D Drag coefficient of an airfoil C L Lift coefficient of an airfoil F Tip-loss factor ω Rotor rotational speed c Blade chord length c Vector of the blade chord lengths α Angle of attack β Pitch control angle β 0 Fixed pitch angle φ Angle of relative wind speed with rotor plane γ Twist angle γ i Actual twist angle for the ith blade element γ i0 Pre-twist angle for the ith blade element γ t Twist angle for the tip of the blade γ r Twist angle for the root of the blade γ γ γ Vector of the blade twist angles σ Solidity ratio k Shape factor of Weibull distribution WT Wind Turbine HAWT Horizontal-Axis Wind Turbine BEMT Blade Element Momentum Theory FPB Fixed Pitch Blade BPC Blade with Pitch Control MB Morphing Blade AEP Annual Energy Production OSU Ohio State University 1 Introduction

Fig. 2 :

 2 Fig. 2: Velocities and forces on a blade element

Fig. 4 :

 4 Fig. 4: C D and C L as a function of α for S809 airfoil

Fig. 5 :Fig. 6 :

 56 Fig. 5: Twist and chord distributions of Phase VI WT blades

Fig. 9 :

 9 Fig. 9: Schematic of a morphing blade with active control

Fig. 10 :Fig. 11 : 3

 10113 Fig. 10: Schematic of a simplified morphing blade

  (27). The lower and upper bounds of the decision variables x are: 0 deg ≤ γ r ≤ 35 deg and -5 deg ≤ γ t ≤ 15 deg.

  mulated as follows: minimize -P(x, p) over x = [γ r γ t ] T p = [b c r t r r ω ρ v re ] subject to P ≤ 19.8 kW 0 deg ≤ γ r ≤ 35 deg -5 deg ≤ γ t ≤ 15 deg (28) 5.2 Optimization Problem for the Blade With Pitch Control For the BPC, the decision variable is the pitch control angle: β. The distribution (γ γ γ) of the pre-twist of the blade is shown in Fig. 5 (shifted by -4.775 deg). The actual twist of the blade is calculated with Eqn. (26). The decision variable β is bounded between -5 deg and 25 deg.

Fig. 12 :

 12 Fig. 12: Optimal performance of the MB and BPC at different wind speeds

Fig. 13 :

 13 Fig. 13: Twist angle as a function of the blade radius for v = 5 m/s, v = 10 m/s and v = 15 m/s

Fig. 15 :

 15 Fig. 15: AEP of the wind turbines as a function of the average wind speed

Table 2 :

 2 Maximum power produced by wind turbines containing FPB, MB and BPC for different wind speeds

		FPB		MB		BPC
	v re	Power Power	γ r	γ t	Power	β
	[m/s]	[kW ]	[kW ]	[deg] [deg]	[kW ]	[deg]
	5.0	2.08	3.28	2.24 0.302	3.36	-0.296
	6.0	3.81	4.96	8.61 -0.002	4.75	1.88
	7.0	5.82	6.41	15.6	-1.78	6.28	2.5
	8.0	7.44	8.1	18.6	-1.28	7.99	4.95
	9.0	8.83	9.74	20.8 0.933	9.69	6.47
	10.0	9.41	11.4	19.9	3.49	11.3	7.85
	11.0	9.48	13.1	27.0	3.24	12.9	9.39
	12.0	9.52	14.9	29.7	4.68	14.4	11.8
	13.0	8.84	16.8	32.2	5.9	15.9	13.7
	14.0	8.09	18.8	33.9	6.51	17.4	15.9
	15.0	8.1	19.8	34.8	1.07	19.0	17.6
	16.0	7.99	19.8	21.2	15.0	19.8	15.4
	17.0	8.12	19.8	21.3	9.26	19.8	13.6
	18.0	8.99	19.8	16.5	12.2	19.8	13.2
	19.0	9.89	19.8	11.5	14.1	19.8	12.8
	20.0	10.8	19.8	19.8	10.2	19.8	12.8
	21.0	11.8	19.8	28.1	8.93	19.8	12.6
	22.0	12.8	19.8	29.4	7.16	19.8	12.0
	23.0	13.8	19.8	16.0	10.0	19.8	10.4
	24.0	14.9	19.8	10.1	11.4	19.8	8.52
	25.0	16.0	19.8	13.1	6.15	19.8	7.09

Table 3 :

 3 AEP of the wind turbines as a function of the average wind speed

		FPB		MB	BPC
	v	AEP	AEP	percentage	AEP	percentage
	[m/s] [MW h/y] [MW h/y]	increase	[MW h/y]	increase
	5.0	24.0	29.9	24.5	29.5	22.7%
	6.0	33.7	43.6	29.4	42.8	27.1%
	7.0	41.6	57.2	37.4	56.0	34.5%
	8.0	48.0	70.1	45.9	68.6	42.7%
	9.0	53.2	81.7	53.5	79.9	50.1%
	10.0	57.5	91.6	59.4	89.6	55.9%
	11.0	60.7	99.4	63.6	97.3	60.2%
	12.0	63.1	105.0	66.4	103.0	63.1%
	13.0	64.5	108.0	68.1	106.0	65.0%
	14.0	65.1	110.0	69.2	108.0	66.2%
	15.0	64.9	110.0	69.7	108.0	66.9%