simulate free fall project, ode45, using reynolds no to get drag coefficient, recursive problem
11 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
LEI Li
el 11 de Sept. de 2021
Comentada: LEI Li
el 12 de Sept. de 2021
It seems this is a recursive problem; I need velocity to calculate Reynolds number;
ODE45 seems no need to use recursive;
now I can't get any result from this code. please give me some advice; thank you.
close all;
clear all;
clc;
tr=[0 15]; %seconds
initv=[0 0]; %start 600 m high
[t,y,Re]=ode45(@rkfalling, tr, initv)
plot(t,y(:,1))
ylabel('x (m)')
xlabel('time(s)')
figure
plot(t,y(:,2))
ylabel('velocity (m/s)')
xlabel('time(s)')
% figure
% plot(t,Re)
% figure
% plot(t,Drag_force)
function [dwdt, Re, Drag_force]=rkfalling(t,w)
g=9.81;
% air density @ standard sea-level condition kg/m3
rou_air = 1.294;
% air viscosity, Pa*s
mu = 17.2e-6;
radius = 0.01;
rou = 7750;
volume_sphere = 4/3*pi*radius^3;
mass =rou*volume_sphere;
%displacement
y = w(1);
%velocity
ydot = w(2);
%Re = w(3);
dwdt=zeros(size(w));
dwdt(1) = ydot;
% drag force
% Reynolds number
Re = rou_air*radius*2*ydot/mu;
%Re=5000;
% drag coefficient
C_D = 24/Re + 2.6*(Re/5.0)/(1+(Re/5.0)^1.52)+0.411*(Re/2.63e5)^(-7.94)/(1+(Re/2.63e5)^(-8.0))+0.25*(Re/1e6)/(1+(Re/1e6));
Drag_force = C_D*0.5*rou_air*pi*radius*2*ydot^2;
%dwdt(2)= -0.2*ydot^2+g;
dwdt(2)= -Drag_force/mass+g;
end
4 comentarios
Respuesta aceptada
Star Strider
el 11 de Sept. de 2021
The NaN values are the result of 0/0 operations, and when that occurs in the first integration, it propagates through all of them.
tr=[0 15]; %seconds
initv=[0 0]+eps; %start 600 m high
[t,y,Re]=ode45(@rkfalling, tr, initv)
plot(t,y(:,1))
ylabel('x (m)')
xlabel('time(s)')
figure
plot(t,y(:,2))
ylabel('velocity (m/s)')
xlabel('time(s)')
% figure
% plot(t,Re)
% figure
% plot(t,Drag_force)
function [dwdt, Re, Drag_force]=rkfalling(t,w)
g=9.81;
% air density @ standard sea-level condition kg/m3
rou_air = 1.294;
% air viscosity, Pa*s
mu = 17.2e-6;
radius = 0.01;
rou = 7750;
volume_sphere = 4/3*pi*radius^3;
mass =rou*volume_sphere;
%displacement
y = w(1);
%velocity
ydot = w(2);
%Re = w(3);
dwdt=zeros(size(w));
dwdt(1) = ydot;
% drag force
% Reynolds number
Re = rou_air*radius*2*ydot/mu;
%Re=5000;
% drag coefficient
C_D = 24/Re + 2.6*(Re/5.0)/(1+(Re/5.0)^1.52)+0.411*(Re/2.63e5)^(-7.94)/(1+(Re/2.63e5)^(-8.0))+0.25*(Re/1e6)/(1+(Re/1e6));
Drag_force = C_D*0.5*rou_air*pi*radius*2*ydot^2;
%dwdt(2)= -0.2*ydot^2+g;
dwdt(2)= -Drag_force/mass+g;
end
.
2 comentarios
Star Strider
el 11 de Sept. de 2021
You wiill need to solve that, since I do not understand what you are doing to the extent that I can correct the code. (My areas of expertise do not include fluid dynamics, unfortunately for this problem.)
Más respuestas (1)
Paul
el 11 de Sept. de 2021
Check your intiial conditions. At t = 0, ydot = 0, which yields Re == 0, which then yields isnan(C_D) == true, and subsequently Draf_force and dwdt(2) are both NaN.
Also, that third output from ode45 isn't what you think it is. It's supposed to be the time of event detections, though no events are explicitly defined for this problem so I'm not really sure what's being returned. But it's certainly not the Reynolds number at each time step. In fact, I'm not even sure that the solver knows or cares about the second and third outputs from rkfalling().
2 comentarios
Paul
el 11 de Sept. de 2021
Once you correct the NaN issue, you can get the Re vs t by calling rkfalling() after you get the solution from ode45
[t,y] = ode45(@rkfalling,tr,initv);
[~,Re,Drag] = rkfalling(t,y);
Ver también
Productos
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!