3. 弹簧-质量系统的自由振动分析

弹簧-质量系统的自由振动分析

本节以单自由度弹簧-质量系统为对象，完成从建模、固有频率与振型计算，到自由振动响应仿真验证的完整流程。系统由质量块 m、线性弹簧 k 和固定支座组成，质量块仅沿水平方向运动，忽略摩擦与阻尼。

取质量块水平位移 x 为广义坐标，原点设在弹簧自然长度处。由牛顿第二定律或拉格朗日方程可得运动方程：

m·ẍ + k·x = 0
令 ωn = √(k/m)  →  ẍ + ωn²·x = 0

该方程的通解为 x(t) = A·cos(ωn·t) + B·sin(ωn·t)，其中 A = x₀，B = v₀/ωn。固有频率 fn = ωn/(2π) 只由 m 和 k 决定，与初始条件无关；振幅和相位则由初始位移 x₀ 与初始速度 v₀ 决定。单自由度系统只有一个振型，即质量块沿弹簧轴线做纯平动，振型向量归一化后为 [1]。

下面用 Python 对上述模型进行数值仿真，并与解析解对比。采用四阶龙格-库塔法积分状态方程 [x, v]' = [v, -ωn²·x]。

import numpy as np
import matplotlib.pyplot as plt

# 系统参数
m = 2.0          # kg
k = 200.0        # N/m
x0 = 0.05        # m
v0 = 0.0         # m/s

wn = np.sqrt(k / m)          # 固有频率 rad/s
fn = wn / (2 * np.pi)        # Hz
T = 1.0 / fn                 # 周期 s

print(f"ωn = {wn:.4f} rad/s")
print(f"fn = {fn:.4f} Hz")
print(f"T  = {T:.4f} s")

# 状态方程: dx/dt = v, dv/dt = -wn^2 * x
def deriv(state):
    x, v = state
    return np.array([v, -wn**2 * x])

# RK4 积分
dt = 0.001
t_end = 2.0
n_steps = int(t_end / dt)
t = np.linspace(0, t_end, n_steps + 1)
state = np.array([x0, v0])
x_num = np.zeros(n_steps + 1)
x_num[0] = state[0]

for i in range(n_steps):
    k1 = deriv(state)
    k2 = deriv(state + 0.5 * dt * k1)
    k3 = deriv(state + 0.5 * dt * k2)
    k4 = deriv(state + dt * k3)
    state = state + (dt / 6.0) * (k1 + 2*k2 + 2*k3 + k4)
    x_num[i + 1] = state[0]

# 解析解
x_exact = x0 * np.cos(wn * t) + (v0 / wn) * np.sin(wn * t)

# 误差
err = np.max(np.abs(x_num - x_exact))
print(f"最大绝对误差 = {err:.2e} m")

# 绘图
plt.figure(figsize=(10, 4))
plt.plot(t, x_exact, 'r-', linewidth=2, label='解析解')
plt.plot(t, x_num, 'b--', linewidth=1.5, label='RK4 数值解')
plt.xlabel('时间 t (s)')
plt.ylabel('位移 x (m)')
plt.title('弹簧-质量系统自由振动响应')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

运行上述代码后，预期输出为：ωn = 10.0000 rad/s、fn = 1.5915 Hz、T = 0.6283 s，数值解与解析解的最大绝对误差应在 1e-6 m 量级（取决于步长）。曲线应呈现等幅余弦振荡，周期约 0.628 s，振幅 0.05 m。