### Junhyuk (view profile)

on 2 Jun 2019
Latest activity Commented on by Junhyuk

on 3 Jun 2019

### Geoff Hayes (view profile)

I need value of 'Wind power energy AEP'
function [P1,P2,P3,P4]=V90
% syntax: function [P1,P2,P3,P4]=V90
% Input of all required parameters of the Vestas V90 wind turbine
% Outputs
% aerodynamic parameters P1=[rho,kp]
% turbine parameters P2=[R,Nb,Jb,kb,mt,dt,kt,nu,Jr,dr,kr,Jg,omg0,kg1,kg2]
% nominal values P4=[Vn,lambdan,thetan]
% air density [kg/m3]
rho=1.25;
% power loss factor [-]; correction factor for the simplifications in BEM [-], see listing 'bem.m'
kp=0.9;
% aerodynamic parameters
P1=[rho,kp];
R=45;
Nb=3;
mb=9600;
% inertia blade (with respect to flapping hinge) [kg m^2]
Jb=3.9e6;
% the stiffness will be determined from the blade flap natural frequency omb [rad/s]
omb=5.71;
kb=Jb*omb^2;
% equivalent mass tower (1/4 mass tower + tower head mass) [kg]
mt=160000;
% stifness tower [N/m]
% the stiffness will be determined from the tower natural frequency [rad/s]
omt=2.51;
kt=mt*omt^2;
% damping tower [N/(m/s)]; 1% critical damping assumed
dt=2*0.01*sqrt(kt*mt);
% inertia generator [kg m^2]
Jg=60;
% nominal terminal (line-to-line) voltage generator [V]
Un=960;
% nominal generator power [W]
Pn=3e6;
% nominal generator shaft angular velocity [rad/s]
omgn=2*pi*25;
% field current generator [A]
If=80;
% number of pole pairs [-]
p=2;
% transmission ratio [-]
nu=87;
% inertia rotor [kg m^2]
Jr=Nb*Jb;
kr=1.6e8;
% total inertia transmission [kg m^2]
Jtot=(nu^2*Jg*Jr)/(nu^2*Jg+Jr);
% damping transmission [Nm/(rad/s)]; 3% critical damping assumed
dr=2*0.03*sqrt(kr*Jtot);
% turbine parameters
P2=[R,Nb,Jb,kb,mt,dt,kt,nu,Jr,dr,kr,Jg,Un,Pn,omgn,If,p];
% number of blade elements [-]
Ns=6;
% radial positions (with respect to rotor axis) of blade sections [m]; not necessary equidistant
% Note: the borders of the blade sections should be given; i.e. Ns+1 values
% first value is start of aerodynamic aerofoil; last value is blade tip (r=R)
r=[4 6.6 10.6 18.5 30.4 41 45];
% chord of blade sections [m]
c=[3.1 3.9 3.9 3.1 2.1 1.3 0.03];
% twist of blade sections [degrees];
% by definition, the last value equals 0 (blade tip)
thetat=[13 13 11 7.8 3.3 0.3 0];
% check
if (length(r) ~= Ns+1) error('number of radial positions not correct');end
if (length(c) ~= Ns+1) error('number of chord values not correct');end
if (length(thetat) ~= Ns+1) error('number of twist values not correct');end
P3=[r;c;thetat];
% nominal wind speed [m/s]
Vn=12;
% nominal tip speed ratio [-]
lambdan=7.8;
% nominal blade pitch angle [degrees]
thetan=-1.5;
% nominal values
P4=[Vn,lambdan,thetan];
--------------------------------------------------------------------------------
function [Dax,Mbeta,Mr,P,Cdax,Cp,a,theta,omr]=powercurve1(windturbine,V)
% syntax: function [Dax,Mbeta,Mr,P,Cdax,Cp,a,theta,omr]=powercurve1(windturbine,V)
% Determination of the characteristics of a VARIABLE SPEED REGULATED wind turbine
% axial force versus wind speed Dax - V
% aerodynamic flap moment versus wind speed Mbeta - V
% aerodynamic rotor torque versus wind speed Mr - V
% aerodynamic power versus wind speed P - V
% thrust coefficient versus wind speed Cdax - V
% power coefficient versus wind speed Cp - V
% induction factor versus wind speed a - V
% blade pitch angle versus wind speed theta - V
% rotor angular velocity versus wind speed omr - V
% It is assumed that the wind turbine has an optimal lambda control, so:
% Partial load (V<=Vn): lambda=lambdan, theta=thetan
% Full load (V>Vn): omr=omrn; theta such that power equals nominal power
%
% Outputs:
% Dax: axial force [N]
% Mbeta: aerodynamic flap moment [Nm]
% Mr: aerodynamic rotor torque [Nm]
% P: aerodynamic power [W]
% Cdax: thrust coefficient [-]
% Cp: power coefficient [-]
% a: induction factor [-]
% theta: blade pitch angle [degrees]
% omr: rotor angular velocity [rad/s]
% Inputs:
% windturbine: name of file with wind turbine parameters (string)
% e.g.: 'LW50'
% V: vector with wind speeds [m/s]
% required parameters
[P1,P2,P3,P4]=feval(windturbine);
R=P2(1);
% transmission ratio [-]
nu=P2(8);
% nominal wind speed [m/s]
Vn=P4(1);
% nominal tip speed ratio [-]
lambdan=P4(2);
% nominal blade pitch angle [degrees]
thetan=P4(3);
% stationary conditions: flap speed and tower top speed equal zero
xd=0;
% nominal rotor angular velocity
omrn=lambdan*Vn/R;
% nominal (mechanical) generator angular velocity
omgn=nu*omrn;
% nominal mechanical power Pn (wind speed equal to nominal wind speed; blade pitch angle equal to
% nominal blade pitch angle; rotor angular velocity equal to nominal rotor angular velocity)
N=length(V);
% calculation of aerodynamic forces, moments etc. for each wind speed
for i=1:N
if V(i) <= Vn
% partial load conditions (wind speed smaller or equal nominal wind speed)
% the tip speed ratio equals nominal tip speed ratio, so the rotor angular velocity equals:
omr(i)=lambdan*V(i)/R;
theta(i)=thetan;
% calculation of the aerodynamic forces, moments etc. by means of blade element-momentum method (BEM)
else
% ful load conditions (wind speed larger than nominal wind speed)
% the rotor angular velocity is kept constant at nominal value
omr(i)=lambdan*Vn/R;
% the blade pitch angle should be such that the power equals nominal power;
% it is assumed that the blade pitch control is to zero-lift
% Use is made of the standard Matlab routine 'fzero' to find a zero of the function 'fun_power.m'; 'fzero' varies
% the blade pitch angle (in the range thetan to 50) until 'fun_power' equals (about) zero.
warning off
options=optimset('Display','off');
theta(i)=fzero('fun_power',[thetan 50],options,V(i),Pn,P1,P2,P3,P4);
warning on
% since the blade pitch angle is determined, the aerodynamic forces, moments and power
% can be calculated by means of the blade element - momentum method (BEM)
end
end
--------------------------------------------------------------------------------
function demo_windsim
K = menu('Choose a demo','cp-lambda','power curve','BEM','cp-lambda (Matlab 6.0 version)','power curve (Matlab 6.0 version)','BEM (Matlab 6.0 version)');
if K==1
playshow cplambda_demo
elseif K==2
playshow powercurve_demo
elseif K==3
playshow bem_demo
elseif K==4
playshow cplambda_demo60
elseif K==5
playshow powercurve_demo60
else
playshow bem_demo60
end
--------------------------------------------------------------------------------
V=1:1:25;
Vavg=5.1;
f=pi/2*V/Vavg^2.*exp(-pi/4*(V/Vavg).^2);
[Dax,Mbeta,Mr,P,Cdax,Cp,a,theta,omr]=powercurve1(V90,V);
figure(1);
plot(V,P)
=======================================================================
But I can't get value of AEP.
this is error code.
??? Error using ==> feval
Argument must contain a string or function_handle.
Error in ==> powercurve1 at 33
[P1,P2,P3,P4]=feval(windturbine);
Error in ==> Untitled at 4
[Dax,Mbeta,Mr,P,Cdax,Cp,a,theta,omr]=powercurve1(V90,V);
How can i get it..?
Thanks for watching my Q. Have a nice day.

### Geoff Hayes (view profile)

on 2 Jun 2019

Junhyuk - try prefixing the V90 (function) parameter with a @. So change the line
[Dax,Mbeta,Mr,P,Cdax,Cp,a,theta,omr]=powercurve1(V90,V);
to
[Dax,Mbeta,Mr,P,Cdax,Cp,a,theta,omr]=powercurve1(@V90,V);
That got me past the same error message. The next error had to do with the bem function being undefined.

Junhyuk

on 3 Jun 2019