MATLAB Answers

simulate free fall project, ode45, using reynolds no to get drag coefficient, recursive problem

조회 수: 11(최근 30일)
LEI Li
LEI Li 2021년 9월 11일
댓글: LEI Li 2021년 9월 12일
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

댓글을 달려면 로그인하십시오.

채택된 답변

Star Strider
Star Strider 2021년 9월 11일
The NaN values are the result of 0/0 operations, and when that occurs in the first integration, it propagates through all of them.
Add eps (or some other small number) to ‘initv’ and it works.
tr=[0 15]; %seconds
initv=[0 0]+eps; %start 600 m high
[t,y,Re]=ode45(@rkfalling, tr, initv)
t = 121×1
1.0e+00 * 0 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.0001 0.0001 0.0003
y = 121×2
0.0000 0.0000 0.0000 0.0001 0.0000 0.0001 0.0000 0.0002 0.0000 0.0002 0.0000 0.0005 0.0000 0.0007 0.0000 0.0010 0.0000 0.0012 0.0000 0.0025
Re = 6×1
30 1 187 0 0 0
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
Star Strider
Star Strider 2021년 9월 11일
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.)

댓글을 달려면 로그인하십시오.

추가 답변(1개)

Paul
Paul 2021년 9월 11일
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

댓글을 달려면 로그인하십시오.

Community Treasure Hunt

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

Start Hunting!

Translated by