3. 简谐激励下的稳态响应解析求解

上一节我们建立了单自由度受迫振动的动力方程，这一节把外力具体化为简谐形式，求出它的稳态解析解。设系统质量 m、阻尼系数 c、刚度 k，外力为 F(t) = F₀ sin(ωt)，则运动方程为

mẍ + cẋ + kx = F₀ sin(ωt)

先统一符号约定：位移 x 以静平衡位置为原点、向下为正，速度、加速度同向取正；外力相位以 sin 为基准。后文所有相位差都按“响应滞后于激励”的角度给出，正值表示滞后。

从通解到稳态解

方程的全解由齐次解与特解叠加而成。齐次解对应自由振动，因阻尼存在会按 e^−ζω_nt 衰减，经过一段时间后消失，称为瞬态分量。特解则与激励同频率持续存在，称为稳态响应。工程上更关心长期运行状态，因此我们只保留稳态项，设

x_p(t) = X sin(ωt − φ)

代入方程，利用 sin 与 cos 的展开并令同频项系数相等，可得幅值与相位：

X = F₀ / √[(k − mω²)² + (cω)²]

tan φ = cω / (k − mω²)，且 0 ≤ φ ≤ π

引入固有频率 ω_n = √(k/m)、阻尼比 ζ = c/(2√(km))、频率比 r = ω/ω_n，并记静位移 X₀ = F₀/k，可写成无量纲形式：

X/X₀ = 1 / √[(1 − r²)² + (2ζr)²]

tan φ = 2ζr / (1 − r²)

这就是幅频特性与相频特性的解析表达式，也是后续数值验证的对照基准。

共振条件与品质因数

观察幅频表达式，分母最小时幅值最大。对 r 求极值可得位移共振频率 r_res = √(1 − 2ζ²)，对应峰值 X_max/X₀ = 1/(2ζ√(1 − ζ²))。当 ζ 很小时，共振近似发生在 r ≈ 1，峰值近似为 1/(2ζ)。这说明阻尼越小，共振峰越高越尖锐，工程上必须避开这一频段。

品质因数 Q = 1/(2ζ) 是衡量共振尖锐程度的常用指标。它还有两个等价含义：一是共振时位移放大倍数约等于 Q；二是半功率带宽 Δr ≈ 1/Q，即幅值降到峰值的 1/√2 时对应的频率区间宽度。例如 ζ = 0.05 时 Q = 10，共振放大十倍，带宽仅约 0.1，系统对频率非常敏感；ζ = 0.5 时 Q = 1，响应平缓，几乎没有明显共振峰。

相位变化规律

相位 φ 随频率比的变化很有规律：r → 0 时 φ → 0，响应与激励几乎同相；r = 1 时无论阻尼大小，φ 恒等于 90°，响应滞后四分之一周期，此时速度与激励同相，输入功率最大；r → ∞ 时 φ → 180°，响应与激励反相。阻尼只影响过渡的快慢：ζ 越小，φ 在 r = 1 附近翻转越陡峭。这一特征在实验中被用来精确识别固有频率——测量相位达到 90° 的频率即为 ω_n，比找幅值峰值更稳健。

数值计算验证

下面用一段 Python 代码把解析式与数值积分结果对照。取 m = 2 kg、k = 800 N/m、ζ = 0.05、F₀ = 10 N，则 ω_n = 20 rad/s，c = 2ζ√(km) = 2 N·s/m，X₀ = 0.0125 m。选激励频率 ω = 18 rad/s（r = 0.9），用四阶龙格库塔积分足够长时间后取稳态段，提取幅值与相位，与解析值比较。

import numpy as np
from scipy.integrate import solve_ivp

m, k, zeta, F0 = 2.0, 800.0, 0.05, 10.0
wn = np.sqrt(k/m); c = 2*zeta*np.sqrt(k*m)
w = 18.0; X0 = F0/k

def ode(t, y):
    x, v = y
    return [v, (F0*np.sin(w*t) - c*v - k*x)/m]

T = 2*np.pi/w
sol = solve_ivp(ode, [0, 60*T], [0, 0], t_eval=np.linspace(0, 60*T, 20000))
t, x = sol.t, sol.y[0]
mask = t > 50*T                      # 只取稳态段
ts, xs = t[mask], x[mask]
amp_num = (xs.max() - xs.min())/2    # 数值幅值
# 解析值
r = w/wn
X_ana = X0/np.sqrt((1-r**2)**2 + (2*zeta*r)**2)
phi_ana = np.arctan2(2*zeta*r, 1-r**2)
print(f"解析幅值 {X_ana:.5f} m, 数值幅值 {amp_num:.5f} m")
print(f"解析相位 {np.degrees(phi_ana):.2f}°")

按上述参数，解析幅值约为 0.0125/√[(1−0.81)² + (2×0.05×0.9)²] ≈ 0.0125/0.206 ≈ 0.0607 m，相位 φ = arctan2(0.09, 0.19) ≈ 25.3°。数值积分得到的稳态幅值应与之一致（差异来自采样离散与瞬态残留，通常在千分之一量级）。若把 ω 改为 20 rad/s，r = 1，解析幅值变为 X₀/(2ζ) = 0.125 m，相位恰为 90°，这是检验代码与公式是否对齐的最好算例。

把解析式、曲线与数值结果三者对齐，就完成了本节的任务闭环：公式给出任意频率下的幅值与相位，曲线展示全局趋势与共振特征，数值计算则验证公式在具体参数下确实成立。后续做仿真时，这套解析结果将作为校核仿真模型是否正确的最直接基准。