代码:
clc
clear
% Define the time span
tspan = [0 10]; % Time span for the simulation (from 0 to 10 seconds)
v0 = 0;
% Solve the differential equation using ode45
[t, v] = ode45(@ode_function, tspan, v0);
% Plot the velocity as a function of time
plot(t, v)
xlabel('Time (s)');
ylabel('Velocity (m/s)');
title('Velocity vs. Time');
grid on;
function dvdt = ode_function(t, v)
% Define parameters
rho_B = 1000; % Density of the body (kg/m^3)
V_B = 0.01; % Volume of the body (m^3)
g = 9.81; % Gravitational acceleration (m/s^2)
rho_W = 1.225; % Density of the fluid (kg/m^3)
A_B = 0.1; % Cross-sectional area of the body (m^2)
m_B = rho_B * V_B; % Mass of the body (kg)
D_B = 0.2; % Diameter of the body (m)
visc_W = 1.789e-5; % Dynamic viscosity of the fluid (N*s/m^2)
% Calculate Reynolds number (Re)
Re = (rho_W * v * D_B) / visc_W;
% Calculate drag coefficient (C_D) using the given formula
C_D = (24/Re) + (2.6 * (Re/5) / (1 + (Re/5)^1.52)) + (0.41 * (Re/263000)^-7.94 / (1 + (Re/263000)^-8)) + (Re^0.8 / 461000);
% Calculate forces
F_B = rho_B * V_B * g;
F_D = (0.5 * rho_W * A_B * C_D) * v^2;
F_W = m_B * g;
% Calculate acceleration (dv/dt)
dvdt = (((F_B - F_D - F_W))/m_B);
end
这段代码应该在使用时输出正确的图,但所有值v
(除了第一个 0 之外)都是未定义/NaN。问题可能出在ode_function()
函数内部,但我对 ode45 不太有经验。
我尝试编辑 dvdt 的定义方式,并在将其更改为更简单的方程时得到了正确的结果,例如,2t
但这是我能够获得有用的结果。
通过使用断点进行调试(您应该在 Matlab 中自行调试)我发现了以下内容:
v0=0
您在 ode_function 中初始化并使用它作为 的因子Re = (rho_W * v * D_B) / visc_W
,因此Re=0
。这导致
C_D = (24/Re) + (2.6 * (Re/5) / (1 + (Re/5)^1.52)) + (0.41 * (Re/263000)^-7.94 / (1 + (Re/263000)^-8)) + (Re^0.8 / 461000)
您除以它的公式(24/Re
-> Inf)。稍后在公式中,您将提升 Re(其中一种情况是Re^0.8
-> Inf)。这会导致Inf/Inf
-> NaN。因此,你之后所做的任何事情都不会产生任何有用的结果,包括情节。
简单的解决方案是不使用 0 作为输入 v0 和/或验证(如 <>0)它作为函数的参数