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.

[Mathieu NOE]

hi

missing function windTurbineRotorModel

[Walter Roberson]

plot(WindSpeeds,generatorTorque(i))

You are after the for i loop at that point, so i will have the scalar value that it was last assigned, which is length(time) . So you are asking to plot() with a vector on the x and a scalar for the y.

[Sebastian Sunny]

I have added my function to the question at the bottom

[Sebastian Sunny]

I thought that was the problem so i removed the (i) on all the plot functions but i end up with a completly different graph.

[Walter Roberson]

if (WindSpeeds < Vcutin) 

Is WindSpeeds a scalar or a vector? If it is a scalar, then that statement is okay. If it is a vector then the statement is questionable.

    elseif all(WindSpeeds >Vcutin) && all(WindSpeeds <Vrated)

all those all() -- is WindSpeeds a vector or a scalar? If it is a scalar then if Vcutin is a vector then you have problems back in the if . If it is a scalar and Vcutin is a scalar then why confuse the reader by using all() ?

Be consistent in how you treat your inputs.

[Sebastian Sunny]

thank you it did end up being a problem with the function ive re written it thank you

[Stephen23]

Original question by Sebastion Sunny retrieved from Google Cache:

https://webcache.googleusercontent.com/search?q=cache:BYIMeK3tlrkJ:https://www.mathworks.com/matlabcentral/answers/1606430-not-sure-whats-wrong-with-the-values-in-this-plot

https://web.archive.org/web/20211210124935/https://www.mathworks.com/matlabcentral/answers/1606430-not-sure-whats-wrong-with-the-values-in-this-plot

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)

Answer 1

[Mathieu NOE]

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