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.
12 次查看(过去 30 天)
显示 更早的评论
% 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
回答(1 个)
darova
2020-4-6
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 个评论
darova
2020-4-6
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);
% ...
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Particle & Nuclear Physics 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!