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

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]

I can’t follow the code, at least in part because my fluid dynamics backgroound is not sufficiently robust.

Some part of ‘w’ or ‘dwdt’ is likely the velocity. Which one is it?

[LEI Li]

w is the vector, represent displacement y and velocity;

in order to get Re, I need velocity;

velocity is coming from function

I can't get the displacement

[darova]

Velocity is ? What is the problem? Don't you have function ?

[LEI Li]

The problem now is I got NAN for y; no errors.

Answer 1

[Star Strider]

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

.

[LEI Li]

The Reynolds number Re is still not right. I thought every time I got a new velocity, or velocity is updated, the Reynolds number should also update.

but here, it is just has six values, and the values seem not right, bc the Re should be greater than 6000;

a steel ball with radius = 5mm falling from sky, terminal velocity should be greater than 4.5m/s; it doesn't make sense.

[Star Strider]

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.)

Answer 2

[Paul]

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().

[Paul]

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);

[LEI Li]

Thank you. I will dig it and find out the solution