Not sure whats wrong with the values in this plot

2 views (last 30 days)
Sebastian Sunny
Sebastian Sunny on 8 Dec 2021
Edited: Stephen23 on 10 Dec 2021
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
Stephen23
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)

Products

Community Treasure Hunt

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

Start Hunting!

Translated by