I can't solve this problem. Please help me

10 visualizzazioni (ultimi 30 giorni)
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.

Risposta accettata

Geoff Hayes
Geoff Hayes il 2 Giu 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.

Più risposte (0)

Categorie

Scopri di più su Multibody Modeling in Help Center e File Exchange

Tag

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by