Not sure whats wrong with the values in this plot

1 view (last 30 days)
Hi guys, ive been stuck on this bit of code for awhile as im getting a empty plots, im not sure of what im doing wrong is to do with the actaul variables being wrong or is i ahve used the for loop vairale wrong.
  7 Comments
Sebastian Sunny
Sebastian Sunny on 9 Dec 2021
thank you it did end up being a problem with the function ive re written it thank you
Stephen23
Stephen23 on 10 Dec 2021
Edited: Stephen23 on 10 Dec 2021
Original question by Sebastion Sunny retrieved from Google Cache:
Not sure whats wrong with the values in this plot
Hi guys, ive been stuck on this bit of code for awhile as im getting a empty plots, im not sure of what im doing wrong is to do with the actaul variables being wrong or is i ahve used the for loop vairale wrong. I' ve attached a picture of what the graph should look like and the result im getting.
Cp = 0.335;
Ct = 0.042;
Vrated = 11.5; %m/s
Vcutin = 3; %m/s
Vcutout = 25; %m/s
D = 171;
k = 6e6;
j = 100e5;
%Create wind,time and rotational speed variables
WindSpeeds = linspace(0,30,30001);%m/s
deltaT = 0.01;
time = 0:deltaT:300;
%preallocation
rotorTorque = zeros(length(WindSpeeds),1); %Nm
turbinePower = zeros(length(WindSpeeds),1);
generatorTorque = zeros(length(WindSpeeds),1);
omegaRotor = zeros(length(WindSpeeds),1);
%eulers method
for i = 2:length(time)
omegaRotor(i) = omegaRotor(i-1) + deltaT*((windTurbineRotorModel(WindSpeeds(i),Ct,D,Vcutout,Vrated,Vcutin))-(k*omegaRotor(i-1).^2)/j);
generatorTorque(i) = (k*(omegaRotor(i).^2));
rotorTorque(i) = windTurbineRotorModel(WindSpeeds(i),Ct,D,Vcutout,Vrated,Vcutin);
end
yyaxis left
plot(WindSpeeds,omegaRotor)
hold on
plot yy
yyaxis right
plot(WindSpeeds,generatorTorque(i))
hold on
yyaxis right
plot(WindSpeeds,rotorTorque)

Sign in to comment.

Accepted Answer

Mathieu NOE
Mathieu NOE on 9 Dec 2021
hello Seb
I revisited your code and found the main bug . In your equations the term that computes the rotor acceleration ( = net torque divided by rotor inertia) was negative and very big - so not physically meaningfull
by splitting the long line
omegaRotor(i) = omegaRotor(i-1) + deltaT*((windTurbineRotorModel(WindSpeeds(i),Ct,D,Vcutout,Vrated,Vcutin))-(k*omegaRotor(i-1).^2)/j);
into smaller bits I saw that the division (by rotor inertia) was applied only on the second term of the net torque , so this was the major hurdle
BTW it took me some time to understand that the "j" in your constant was indeed refering to the rotor inertia, so I prefer to use a more explicit variable name (which is one of the programmer's good practice , beside not using i and j which are reserved for complex numbers)
otherwise a few minor things could be upgraded , I let you go through this version of the code
also I didn't put much effort to label the figure, my apologize, I asssume you can do it by yourself
all the best
clc
clearvars
Cp = 0.335;
Ct = 0.042;
Vrated = 11.5; %m/s
Vcutin = 3; %m/s
Vcutout = 25; %m/s
D = 171;
k = 6e6;
Rotor_Inertia = 100e5; % do not use i or j !!
%Create wind,time and rotational speed variables
WindSpeeds = linspace(0,30,3001);%m/s % found out 301 or 3001 samples is way enough...
time = linspace(0,300,length(WindSpeeds));
deltaT = mean(diff(time));
%preallocation
rotorTorque = zeros(length(WindSpeeds),1); %Nm
turbinePower = zeros(length(WindSpeeds),1);
generatorTorque = zeros(length(WindSpeeds),1);
omegaRotor = zeros(length(WindSpeeds),1);
%eulers method
for ci = 2:length(time)
rotorTorque(ci) = windTurbineRotorModel(WindSpeeds(ci),Ct,D,Vcutout,Vrated,Vcutin);
net_torque = (rotorTorque(ci) - (k*omegaRotor(ci-1).^2));
accel = net_torque/Rotor_Inertia;
omegaRotor(ci) = omegaRotor(ci-1) + deltaT*accel;
generatorTorque(ci) = (k*(omegaRotor(ci).^2));
rotorTorque(ci) = windTurbineRotorModel(WindSpeeds(ci),Ct,D,Vcutout,Vrated,Vcutin);
end
% yyaxis left
plot(WindSpeeds,omegaRotor)
yyaxis right
plot(WindSpeeds,generatorTorque(ci))
yyaxis right
plot(WindSpeeds,rotorTorque)
%%%%%%%%%%%%%%%%%%%%%%
function [rotorTorque] = windTurbineRotorModel(WindSpeeds,Ct,D,Vcutout,Vrated,Vcutin)
if (WindSpeeds <= Vcutin)
rotorTorque = 0;
% elseif all(WindSpeeds >Vcutin) && all(WindSpeeds <Vrated)
elseif (WindSpeeds >Vcutin) && (WindSpeeds <Vrated)
rotorTorque = (Ct)*(1/2)*(1.23)*pi*((D/2)^3)*WindSpeeds^2;
% elseif all(WindSpeeds >Vrated) && all(WindSpeeds < Vcutout)
elseif (WindSpeeds >=Vrated) && (WindSpeeds <= Vcutout)
rotorTorque = (Ct)*(1/2)*(1.23)*pi*((D/2)^3)*Vrated^2;
else
rotorTorque = 0;
end
end

More Answers (0)

Categories

Find more on MATLAB in Help Center and File Exchange

Products

Community Treasure Hunt

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

Start Hunting!