function dhdt = taper_ode(t, h, r_t, r_b, H, A_out)
g = 9.81;
% 计算当前高度 h 的半径 r
r = r_b + (r_t - r_b) * (h / H);
% 计算当前水面面积 A_surface
A_surface = pi * r^2;
% 计算 dh/dt
dhdt = - (A_out * sqrt(2 * g * h)) / A_surface;
% 确保水排空后 dh/dt 为 0
if h <= 0
dhdt = 0;
end
end
%% 水桶排空仿真脚本
% --- 1. 定义参数 ---
% 几何参数必须手动选择,以确保体积 V 约为 1m^3 (1000升)。
% V = (1/3) * pi * H * (r_t^2 + r_b^2 + r_t*r_b)
% A 桶参数 (上宽下窄)
H_A = 1.12; % 初始高度
r_t_A = 0.6; % 顶部半径
r_b_A = 0.1; % 底部半径
% B 桶参数 (上窄下宽)
H_B = 1.12; % 初始高度
r_t_B = 0.1; % 顶部半径
r_b_B = 0.6; % 底部半径
A_out = 0.001; % 水龙头面积 (m^2)
tspan = [0 500]; % 仿真时间 (秒)
h0_A = H_A; % A 桶初始高度
h0_B = H_B; % B 桶初始高度
% --- 2. 求解微分方程 ---
% 求解 A 桶:使用匿名函数将额外的参数传递给 taper_ode
[t_A, h_A] = ode45(@(t, h) taper_ode(t, h, r_t_A, r_b_A, H_A, A_out), tspan, h0_A);
% 求解 B 桶
[t_B, h_B] = ode45(@(t, h) taper_ode(t, h, r_t_B, r_b_B, H_B, A_out), tspan, h0_B);
% --- 3. 分析结果 ---
% 找到排空时间 T (寻找高度 h 第一次小于一个极小值 0.001m 的时间点)
T_A = t_A(find(h_A < 0.001, 1, 'first'));
T_B = t_B(find(h_B < 0.001, 1, 'first'));
disp('---------------------------------');
disp(['A 桶排空时间 (秒): ', num2str(T_A)]);
disp(['B 桶排空时间 (秒): ', num2str(T_B)]);
disp('---------------------------------');
% --- 4. 绘图 ---
figure;
plot(t_A, h_A, 'r', 'LineWidth', 2);
hold on;
plot(t_B, h_B, 'b--', 'LineWidth', 2);
grid on;
legend('桶 A (上宽下窄)', '桶 B (上窄下宽)', 'Location', 'best');
xlabel('时间 (秒)');
ylabel('水面高度 (米)');
title('水桶排空仿真 (体积相同)');
【 在 tsuld 的大作中提到: 】
: 谢谢,方便时可以贴一下仿真代码学习一下。
: 本题确实是A快。逻辑推理如下:
: 假设最终B先流完,由于一开始B下降的液面就比A快那么一丢丢,造成A的出水流速高于B,由于整个过程是连续的,很难想象在B一开始领先,最后仍然领先的过程中会出现A的液面能追上B的情况,也就是全过程A的液面都是高于B的,但这会导致A的出水速度全程大于B,在体积相同时反而后流完,矛盾。故B不可能先流完。
: ...................
--
FROM 175.152.121.*