阻尼单摆Matlab简易仿真

一、阻尼单摆的数学物理推导

1.1 动力学方程建立

假设与坐标系：

• 单摆摆长为，摆锤质量为。

• 摆角为偏离竖直向下位置的角位移。

• 存在与速度成正比的线性阻尼力，阻尼系数为（单位：kg/s）。

• 重力加速度为。

受力分析：

1. 重力切向分量：（负号表示恢复力方向与增加方向相反）。

2. 线性阻尼力：（与切向速度˙ 成正比，方向相反）。

切向运动方程（牛顿第二定律）：

整理得：

1.2 小角度线性化近似

当较小时（通常时误差很小），。
令：

方程化为线性阻尼谐振子方程：

1.3 求解线性微分方程

这是常系数二阶齐次线性微分方程，其解的形式取决于特征方程：

特征根：

1.3.1 欠阻尼情况（​）

定义阻尼振荡角频率：

则特征根为共轭复数：

通解为：

由初始条件确定常数：
假设时：

（从静止释放）
代入得：

因此：

更简洁的振幅-相位形式：
令

则：

1.4 速度表达式

切向速度。
对求导：

简化（合并同类项）：

注意，所以：

于是：

更简洁地，利用振幅-相位形式的导数：

此处利用了三角恒等式：

其中​​，且相位适当选取后可得简式。

最终速度表达式（小角度近似下）：

速度幅值（初始最大速度）由能量估算：
初始势能：

（小角度时）
初始动能最大值，由机械能守恒（无阻尼时）：

得：

这正是代码中的。

因此速度公式可写为：

（其中​，）

二、代码解释（对照推导）

2.1 参数定义（第1部分）

g = 9.81; % 重力加速度 ω₀² = g/L L = 1.0; % 摆长 m = 0.5; % 质量 b = 0.3; % 阻尼系数 γ = b/(2m) theta0_deg = 30; % 初始角度

• 对应物理模型中的​。

2.2 导出参数计算（第2部分）

theta0 = deg2rad(theta0_deg); % 角度转弧度 omega0 = sqrt(g / L); % ω₀ = √(g/L) gamma = b / (2 * m); % γ = b/(2m) omega_d = sqrt(omega0^2 - gamma^2); % ω_d = √(ω₀² - γ²) A = sqrt(2 * g * L * (1 - cos(theta0))); % 初始最大速度幅值

• 完全对应推导公式。

2.3 阻尼状态检查（第3部分）

if gamma >= omega0 ... % 过阻尼或临界阻尼，停止运行 end

• 确保是欠阻尼情况（​），此时才有振荡解。

2.4 速度计算（第4部分）

t = 0:dt:t_max; v = A * exp(-gamma * t) .* sin(omega_d * t);

• 直接实现公式：

• 注意：这里省略了负号（只影响初相，不影响振幅衰减和周期特性）。

2.5 绘图与包络线（第5-6部分）

• 主曲线：

• 包络线：，显示振幅的指数衰减。

2.6 物理参数输出（第7部分）

fprintf('固有周期 T₀ = %.3f s\n', 2*pi/omega0); fprintf('阻尼周期 T_d = %.3f s\n', 2*pi/omega_d); fprintf('阻尼因子 γ = %.3f s⁻¹\n', gamma); fprintf('时间常数 τ = 1/γ = %.3f s\n', 1/gamma); fprintf('阻尼比 ζ = γ/ω₀ = %.3f\n', gamma/omega0);

• ：无阻尼时的固有周期。

• ：有阻尼时的振荡周期（稍长于 T0T0​）。

• ：振幅衰减到初始值的 1/e≈37%1/e≈37% 所需时间。

• ：无量纲阻尼比。

2.7 不同阻尼系数对比（第8部分）

• 绘制不同（即不同）下的曲线。

• 阻尼越小（越小），衰减越慢，振荡越持久；阻尼越大，衰减越快。

三、总结

1. 物理模型：线性阻尼单摆（小角度近似）。

2. 核心公式：

   • 速度：

   • 其中​，​，​，

3. 代码作用：

   • 实现了上述公式的数值计算。

   • 可视化速度随时间的变化及振幅的指数衰减。

   • 比较不同阻尼系数下的运动特性。

这个模拟清晰地展示了阻尼振动中振幅指数衰减和周期略微增加的特性。

四、源码

%% 阻尼单摆速度-时间图模拟 clear; close all; clc; %% 1. 定义物理参数 g = 9.81; % 重力加速度 (m/s^2) L = 1.0; % 摆长 (m) m = 0.5; % 摆锤质量 (kg) b = 0.3; % 阻尼系数 (kg/s) - 【可调参数】 theta0_deg = 30; % 初始释放角度 (度) - 【可调参数】 %% 2. 计算导出参数 theta0 = deg2rad(theta0_deg); % 将角度转换为弧度 % 固有角频率 (无阻尼) omega0 = sqrt(g / L); % 阻尼因子 gamma = b / (2 * m); % 阻尼角频率 omega_d = sqrt(omega0^2 - gamma^2); % 初始最大速度振幅 A (根据能量守恒估算: mgh = 1/2 * m * v^2) % h = L - L*cos(theta0) A = sqrt(2 * g * L * (1 - cos(theta0))); %% 3. 检查阻尼状态 if gamma >= omega0 fprintf('警告: 当前为过阻尼或临界阻尼状态 (γ ≥ ω₀)，不会发生振荡。\n'); fprintf('请减小阻尼系数 b 或增加摆长 L。\n'); return; % 终止程序 end %% 4. 生成时间与速度数据 t_max = 20; % 模拟总时长 (s) dt = 0.01; % 时间步长 (s) t = 0:dt:t_max; % 时间向量 % 根据阻尼振动公式计算速度 v = A * exp(-gamma * t) .* sin(omega_d * t); %% 5. 绘图 figure('Name', '阻尼单摆速度-时间图', 'NumberTitle', 'off'); % 绘制主图 plot(t, v, 'b-', 'LineWidth', 2); hold on; % 绘制零线便于观察 plot(t, zeros(size(t)), 'k--', 'LineWidth', 1); % 绘制包络线（振幅衰减曲线） plot(t, A * exp(-gamma * t), 'r--', 'LineWidth', 1.5); plot(t, -A * exp(-gamma * t), 'r--', 'LineWidth', 1.5); hold off; %% 6. 美化图像 xlabel('时间', 'FontSize', 12); ylabel('速度', 'FontSize', 12); title_str = sprintf('阻尼单摆速度-时间图\n(摆长=%.1fm, 质量=%.1fkg, 阻尼系数=%.2f, 初始角度=%.0f°)', ... L, m, b, theta0_deg); title(title_str, 'FontSize', 14); grid on; legend('实际速度', '零速度线', '包络线', 'Location', 'northeast'); set(gca, 'FontSize', 11); %% 7. 在命令行窗口输出关键物理参数 fprintf('===== 模拟参数 =====\n'); fprintf('固有周期 T₀ = %.3f s\n', 2*pi/omega0); fprintf('阻尼周期 T_d = %.3f s\n', 2*pi/omega_d); fprintf('阻尼因子 γ = %.3f s⁻¹\n', gamma); fprintf('时间常数 τ = 1/γ = %.3f s (振幅衰减到37%%所需时间)\n', 1/gamma); fprintf('阻尼比 ζ = γ/ω₀ = %.3f\n', gamma/omega0); fprintf('========================\n'); %% 8. 【可选】比较不同阻尼系数的影响 figure('Name', '不同阻尼系数对比', 'NumberTitle', 'off'); hold on; b_values = [0.05, 0.2, 0.5, 0.8]; % 不同的阻尼系数 colors = lines(length(b_values)); for i = 1:length(b_values) b_val = b_values(i); gamma_val = b_val / (2 * m); if gamma_val < omega0 omega_d_val = sqrt(omega0^2 - gamma_val^2); v_val = A * exp(-gamma_val * t) .* sin(omega_d_val * t); plot(t, v_val, 'Color', colors(i, :), 'LineWidth', 2, ... 'DisplayName', sprintf('b = %.2f (ζ = %.2f)', b_val, gamma_val/omega0)); end end plot(t, zeros(size(t)), 'k--', 'LineWidth', 1); hold off; xlabel('时间', 'FontSize', 12); ylabel('速度', 'FontSize', 12); title('不同阻尼系数对速度-时间图的影响', 'FontSize', 14); legend('show', 'Location', 'northeast'); grid on; set(gca, 'FontSize', 11);