【弹簧】多自由度模拟弹簧质量系统Matlab实现

matlab科研助手

于 2025-02-11 11:08:11 发布

阅读量1k

点赞数 28

CC 4.0 BY-SA版权

文章标签： matlab 开发语言

本文链接：https://blog.youkuaiyun.com/matlab_dingdang/article/details/145567117

✅作者简介：热爱科研的Matlab仿真开发者，擅长数据处理、建模仿真、程序设计、完整代码获取、论文复现及科研仿真。

🍎个人主页：Matlab科研工作室

🍊个人信条：格物致知。

🔥 内容介绍

本文旨在阐述如何通过数值方法计算和仿真单自由度 (Single Degree of Freedom, SDOF) 和双自由度 (Two Degree of Freedom, 2DOF) 弹簧质量系统的位置和速度。为避免冗余，本文将略去具体的数学公式推导，相关数学背景知识已在配套的演示幻灯片中详细呈现。读者应参照幻灯片以获得完整的数学基础。本文重点在于讨论数值方法的应用，并展示仿真结果，以深入理解这两种系统的动态特性。

1. 引言

弹簧质量系统是机械工程、土木工程以及其他工程领域中常用的理想化模型，用于分析振动和动态响应。单自由度和双自由度系统是理解更复杂系统行为的基础。精确地预测这些系统的运动状态，对于结构设计、控制系统开发和故障诊断至关重要。解析解虽然可以在某些简单情况下获得，但对于更复杂的情况，例如非线性弹簧或阻尼，以及复杂的激励函数，数值方法成为不可或缺的工具。

2. 数值方法选择

多种数值方法可以用于求解弹簧质量系统的运动方程。常用的方法包括：

欧拉法 (Euler Method): 一种简单的一阶方法，易于实现，但精度较低，且对步长敏感。
改进欧拉法 (Improved Euler Method): 也称为 Heun's method，精度高于欧拉法，但计算量稍大。
龙格-库塔法 (Runge-Kutta Methods): 一系列高阶方法，包括经典的四阶龙格-库塔法 (RK4)，具有较高的精度和稳定性，是常用的选择。
Newmark 方法: 针对结构动力学问题开发的隐式积分方法，适用于求解具有刚性特性的系统。

由于 RK4 方法在精度和计算效率之间取得了良好的平衡，本文选择 RK4 方法作为主要的数值求解方法。

3. 单自由度弹簧质量系统的数值仿真

3.1 模型描述

单自由度弹簧质量系统由一个质量块 (m)、一个弹簧 (k) 和一个阻尼器 (c) 组成。质量块在单一方向上运动。其运动方程可以表示为：

m*x''(t) + c*x'(t) + k*x(t) = F(t)

其中，x(t) 是质量块的位置，x'(t) 是速度，x''(t) 是加速度，F(t) 是外力。

3.2 数值实现

为了使用 RK4 方法，需要将二阶微分方程转换为一阶微分方程组。令 v(t) = x'(t)，则有：

x'(t) = v(t)
v'(t) = (F(t) - c*v(t) - k*x(t)) / m

现在，可以使用 RK4 方法迭代计算 x(t) 和 v(t)。RK4 方法的通用形式如下：

k1 = h * f(t_n, y_n)
k2 = h * f(t_n + h/2, y_n + k1/2)
k3 = h * f(t_n + h/2, y_n + k2/2)
k4 = h * f(t_n + h, y_n + k3)
y_{n+1} = y_n + (k1 + 2*k2 + 2*k3 + k4) / 6

其中，h 是时间步长，y_n 是 t_n 时刻的状态变量，f(t, y) 是微分方程的右侧。在本例中，y = [x, v]，f(t, y) = [v, (F(t) - cv - kx) / m]。

3.3 仿真结果与分析

通过选择合适的 m、k、c 值和外力 F(t)，可以模拟各种不同的系统行为。例如，当 F(t) = 0 且 c = 0 时，系统将发生自由振动。改变阻尼系数 c 可以观察欠阻尼、临界阻尼和过阻尼状态。施加正弦激励可以研究系统的共振现象。

仿真结果通常以位置-时间曲线和速度-时间曲线的形式展示。可以进一步分析系统的固有频率、阻尼比、振幅响应等特征。需要注意的是，选择合适的时间步长 h 对于确保数值解的精度和稳定性至关重要。过大的步长可能导致误差累积甚至发散。

4. 双自由度弹簧质量系统的数值仿真

4.1 模型描述

双自由度弹簧质量系统由两个质量块 (m1, m2)、三个弹簧 (k1, k2, k3) 和两个阻尼器 (c1, c2) 组成。两个质量块分别在各自的方向上运动。其运动方程可以表示为：

m1*x1''(t) + c1*x1'(t) + k1*x1(t) + k2*(x1(t) - x2(t)) = F1(t)
m2*x2''(t) + c2*x2'(t) + k3*x2(t) + k2*(x2(t) - x1(t)) = F2(t)

其中，x1(t) 和 x2(t) 分别是两个质量块的位置，F1(t) 和 F2(t) 是作用在两个质量块上的外力。

4.2 数值实现

类似地，将二阶微分方程转换为一阶微分方程组。令 v1(t) = x1'(t) 和 v2(t) = x2'(t)，则有：

x1'(t) = v1(t)
x2'(t) = v2(t)
v1'(t) = (F1(t) - c1*v1(t) - k1*x1(t) - k2*(x1(t) - x2(t))) / m1
v2'(t) = (F2(t) - c2*v2(t) - k3*x2(t) - k2*(x2(t) - x1(t))) / m2

现在，可以使用 RK4 方法迭代计算 x1(t)、x2(t)、v1(t) 和 v2(t)。在本例中，y = [x1, x2, v1, v2]，f(t, y) 包含上述四个一阶微分方程的右侧。

本文阐述了如何使用数值方法（特别是 RK4 方法）计算和仿真单自由度和双自由度弹簧质量系统的位置和速度。通过仿真，可以深入理解这些系统的动态特性，例如固有频率、阻尼比、振动模式和能量传递。这些仿真结果对于结构设计、控制系统开发和故障诊断具有重要的意义。读者应结合演示幻灯片中的数学推导，以更全面地理解本文所描述的数值方法和仿真结果。未来的工作可以扩展到更复杂的系统，例如多自由度系统和非线性系统，以及更先进的数值方法，例如有限元方法。

📣 部分代码

clc;clear all;close all;global m k w_f F0mass = "Mass of body: ";m = input(mass);sprnConst = "Spring constant: ";k = input(sprnConst);for_w = "Angular frequency of exciting force: ";w_f = input(for_w);F_amp = "Amplitude of exciting force: ";F0 = input(F_amp);dt = .005;t = 0:dt:3;c1 = "Initial velocity: ";vel0 = input(c1);c2 = "Iniitial displacement: ";disp0 = input(c2);y0 = [vel0 disp0]; %[vel disp][tsol,ysol] = ode45('odetest2',0:dt:3,y0);plot(t,ysol(:,2));xlabel('Displacement')ylabel('metres')figureplot(t,ysol(:,1));xlabel('Velocity')ylabel('metre/sec')T = 2*pi*sqrt(m/k)figure% creates mass plotmass = plot(ysol(1,2), 0, 'bo', 'MarkerSize', 30, 'MarkerFaceColor', 'b');% sets axis limitsaxis([min(ysol(:,2))-0.1 max(ysol(:,2))+0.1 -0.1 0.1]);% loops through time steps and update position of masstime = t;for i = 1:length(time)    set(mass, 'XData', ysol(i,2));    pause(dt);end