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]
% blade geometry P3=[r;c;thetat]
% 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];% rotor radius [m]
R=45;% number of blades [-]
Nb=3;% blade mass [kg]
mb=9600;% inertia blade (with respect to flapping hinge) [kg m^2]
Jb=3.9e6;% stifness flap spring [Nm/rad]
% 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;% stiffness transmission [Nm/rad]
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');endif (length(c) ~= Ns+1) error('number of chord values not correct');endif (length(thetat) ~= Ns+1) error('number of twist values not correct');end% blade geometry
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);% rotor radius
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
betad=0;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)
[Dax,Mbeta,Mr,Pn,Cdax,Cp,a]=bem(Vn,thetan,betad,omrn,xd,P1,P2,P3);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; % the blade pitch angle equals nominal blade pitch angle
theta(i)=thetan; % calculation of the aerodynamic forces, moments etc. by means of blade element-momentum method (BEM)
[Dax(i),Mbeta(i),Mr(i),P(i),Cdax(i),Cp(i),a(i)]=bem(V(i),theta(i),betad,omr(i),xd,P1,P2,P3); 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)
[Dax(i),Mbeta(i),Mr(i),P(i),Cdax(i),Cp(i),a(i)]=bem(V(i),theta(i),betad,omr(i),xd,P1,P2,P3); endend--------------------------------------------------------------------------------function demo_windsimK = 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_demoelseif K==2 playshow powercurve_demoelseif K==3 playshow bem_demoelseif K==4 playshow cplambda_demo60elseif K==5 playshow powercurve_demo60else playshow bem_demo60end--------------------------------------------------------------------------------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 ==> fevalArgument 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.
Best Answer