MATLAB Answers

Function returns a vector of length 6, but the length of initial conditions vector is 34. The vector returned by ICEMODELFUNC and the initial conditions vector must have the same number of elements.

4 views (last 30 days)
AeroEng
AeroEng on 6 Apr 2020
Commented: darova on 6 Apr 2020
% Flat Terrain
thetaDEG = 0:45:315; % Degrees
thetaRAD = deg2rad(thetaDEG); % Radians
r = 55; % Meters
m = 0.5; % kg
rhoice = 900; % kg/m^3 (hard rime or glaze)
omegaRPM = 14.5; % rpm
omega = 1.5184364475; % rad/s
Cd = 1;
rhoair = 1.25; % kg/m^3
u = 100; % Meters
l = (m/rhoice)^(1/3);
A = l^2; % m^2
% A = 0.0135; % m^2 - Assumed to be a CUBE shape
Z0 = 0.01; % Meters - Roughness Length
Zh = 125; % Meters - Hub Height
U = 10; % m/s - Velocity at U
g = 9.81; % m/s
% Time
tdelta = 0.2; % seconds
t0 = 0; % seconds
tf = 10; % seconds
t = t0:tdelta:tf;
% Turbine Location
x = 0;
y = 0;
z = 0;
pos = [x,y,z];
% Aerodynamic Drag
Fd = -Cd*rhoair*A*U^2;
for i = 1:length(thetaRAD)
x1 = 0; % x0
x2 = 0; % Vx0
x3(i) = r*cos(thetaRAD(i)); % y0
x4(i) = -r*omega*sin(thetaRAD(i)); % Vy0
x5(i) = r*sin(thetaRAD(i)); % z0
x6(i) = r*omega*cos(thetaRAD(i)); % Vz0
end
% ODE Function
initial = [x1,x2,x3,x4,x5,x6];
tspan = t0:tdelta:tf;
[t,x] = ode45(@IceModelFunc, tspan, initial);
plot(x(1),x(3))
grid on
grid minor
---------------------------------------MY FUNCTION IS BELOW---------------------------------------------
function funcxyz = IceModelFunc(t,x)
m = 0.5;
rhoair = 1.25;
Cd = 1;
A = 0.0068;
U = 10;
g = 9.81;
% x = x1;
% xd = x2;
% y = x3;
% yd = x4;
% z = x5;
% zd = x6;
funcxyz = zeros(6,1);
funcxyz(1) = x(2);
funcxyz(2) = -(1/(2*m))*rhoair*Cd*A*(x(2)-U)*sqrt((x(2)-U)^2*x(4)^2*x(6)^2);
funcxyz(3) = x(4);
funcxyz(4) = -(1/(2*m))*rhoair*Cd*A*x(4)*sqrt((x(2)-U)^2*x(4)^2*x(6)^2);
funcxyz(5) = x(6);
funcxyz(6) = -g - (1/(2*m))*rhoair*Cd*A*(x(5))*sqrt((x(2)-U)^2*x(4)^2*x(6)^2);
end

  3 Comments

AeroEng
AeroEng on 6 Apr 2020
Yes, So I am trying to perform multiple iterations in order to plot distance along the x-direction (which is x(1)) vs distance along the y-direction (which is x(3)) all with different values of theta. I have 8 values for theta and I want to plot all 8 of those values. So yes I have multiple iterations. However the output of the function gives me an vector length of 6 (x(1), x(2), ... x(6)) and my for loop gives me an output for 8 theta values but for all 6 (first two being equal to zero so those are a single value) values so I end up with 34 as a length for the vector out of my for loop. I need to somehow make those equal so I can plot this. I should end up with 8 lines in my plot at the end but I dont get anything due to this issue. Let me know if there is any more confusion I will clear it up. THANK YOU SO MUCH!
AeroEng
AeroEng on 6 Apr 2020
Also just to let you know:
x(1) = Position in x-direction - x
x(2) = Velocity in x-direction - x_dot or dx/dt
x(3) = Position in y-direction - y
x(4) = Velocity in y-direction - y_dot or dy/dt
x(5) = Position in z-direction - z
x(6) = Velocity in z-direction - z_dot or dz/dt

Sign in to comment.

Answers (1)

darova
darova on 6 Apr 2020
You want several solution for thetaRAD
So just put your ode45 function inside for loop
tspan = t0:tdelta:tf;
for i = 1:length(thetaRAD)
x1 = 0; % x0
x2 = 0; % Vx0
x3 = r*cos(thetaRAD(i)); % y0
x4 = -r*omega*sin(thetaRAD(i)); % Vy0
x5 = r*sin(thetaRAD(i)); % z0
x6 = r*omega*cos(thetaRAD(i)); % Vz0
% ODE Function
initial = [x1,x2,x3,x4,x5,x6];
[t,x] = ode45(@IceModelFunc, tspan, initial);
figure(i)
title(sprintf('theta: %d',thetaDEG(i)))
line(t,x)
end

  6 Comments

Show 3 older comments
AeroEng
AeroEng on 6 Apr 2020
How do you mean. I only see a few differences that when i fixed them, it didnt change the outcome. I forgot to comment out the x=x1;xd=x2;...etc. Did that and didnt change it. Also i changed the xx from the ode back to a single x. Lastly I changes the zeros(6,6) back to zeros(6,1). None of that fixed it. What am I missing?
darova
darova on 6 Apr 2020
I mean declaration
If you want to use x1 instead of x(1) declare them inside your function
function func1 = myfuncx(t,x1,x2,x3,x4,x5,x6)
m = 0.5;
rhoair = 1.25;
Cd = 1;
A = 0.0068;
U = 10;
x1 = x(1);
x2 = x(2);
% ...

Sign in to comment.

Sign in to answer this question.


Translated by