10 views (last 30 days)

Show older comments

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

Star Strider
on 11 Sep 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

.

Star Strider
on 11 Sep 2021

Paul
on 11 Sep 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().

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

Start Hunting!