I can't solve this problem. Please help me

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');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

% 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);
  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.

Answer 1

[Geoff Hayes]

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]

Thank you so much! Your answer was very very helpful.