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

Question

LEI Li 2021-9-11

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1451014-simulate-free-fall-project-ode45-using-reynolds-no-to-get-drag-coefficient-recursive-problem

评论： LEI Li 2021-9-12

采纳的回答： Star Strider

在 MATLAB Online 中打开

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 个评论
显示 2更早的评论隐藏 2更早的评论

darova 2021-9-11

Velocity is

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

?

LEI Li 2021-9-11

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

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2021-9-11

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1451014-simulate-free-fall-project-ode45-using-reynolds-no-to-get-drag-coefficient-recursive-problem#answer_785494

在 MATLAB Online 中打开

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 个评论
显示无隐藏无

LEI Li 2021-9-11

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

请先登录，再进行评论。

Answer 2

Paul 2021-9-11

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1451014-simulate-free-fall-project-ode45-using-reynolds-no-to-get-drag-coefficient-recursive-problem#answer_785489

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 个评论
显示无隐藏无

Paul 2021-9-11

在 MATLAB Online 中打开

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 2021-9-12

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

请先登录，再进行评论。

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

4 个评论
显示 2更早的评论隐藏 2更早的评论

采纳的回答

2 个评论
显示无隐藏无

更多回答（1 个）

2 个评论
显示无隐藏无

另请参阅

类别

标签

产品

Community Treasure Hunt

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

4 个评论 显示 2更早的评论隐藏 2更早的评论

采纳的回答

2 个评论 显示 无隐藏 无

更多回答（1 个）

2 个评论 显示 无隐藏 无

另请参阅

类别

标签

产品

Community Treasure Hunt

4 个评论
显示 2更早的评论隐藏 2更早的评论

2 个评论
显示无隐藏无

2 个评论
显示无隐藏无