4. 一维振动系统的数值模拟

弹簧-质量-阻尼系统是振动数值模拟的“最小可用模型”，它把连续体的复杂性剥离，只留下惯性、弹性与耗散三种基本效应。单自由度系统的运动方程为

m·ü(t) + c·u̇(t) + k·u(t) = F(t)

其中 m 为质量，c 为阻尼系数，k 为刚度，F(t) 为外激励。约定位移 u 以静平衡位置为原点、向上为正，速度与加速度同向取正，全文与图示保持一致。数值模拟的任务，是把这条连续时间的微分方程离散成时间步序列，逐步推进求解 u_n。本节围绕两种离散策略——显式积分与隐式积分——展开，分析它们的稳定性与精度，并用自由振动与受迫振动两个案例说明如何选型与判读结果。

显式积分：中心差分法的推进逻辑

显式方法最常用的是中心差分（Central Difference Method）。它用相邻三步的位移近似加速度与速度：

ü_n ≈ (u_{n+1} - 2u_n + u_{n-1}) / Δt²
u̇_n ≈ (u_{n+1} - u_{n-1}) / (2Δt)

代入运动方程并整理，可直接解出下一步位移：

u_{n+1} = [ (2m - kΔt²)·u_n - (m - cΔt/2)·u_{n-1} + F_n·Δt² ] / (m + cΔt/2)

注意右端只含已知量 u_n、u_n−1 与当前荷载 F_n，无需解方程组，每步计算量极小，这是显式方法在大型有限元程序中被广泛采用的根本原因。代价是条件稳定性：中心差分要求时间步长满足

Δt ≤ Δt_cr = 2 / ω_max · ( √(1+ξ²) - ξ )  ≈ 2/ω_max

其中 ω_max 是系统最高固有圆频率，ξ 为阻尼比。对单自由度系统 ω_max = √(k/m)。举例：取 m = 1 kg、k = 100 N/m，则 ω = 10 rad/s，临界步长约 0.2 s；若取 Δt = 0.25 s 超过临界值，数值解会指数发散，即使物理系统本身是衰减的。这一现象必须通过试算识别：位移幅值逐周期放大而非衰减，就是步长越界的典型信号。

隐式积分：Newmark 法的无条件稳定

隐式方法以 Newmark-β 家族为代表，它假设步内加速度线性变化，用两个参数 β、γ 控制精度与稳定性：

u_{n+1} = u_n + Δt·u̇_n + Δt²[(1/2 - β)·ü_n + β·ü_{n+1}]
u̇_{n+1} = u̇_n + Δt[(1 - γ)·ü_n + γ·ü_{n+1}]

与运动方程联立后，得到关于 u_n+1 的线性方程，其等效刚度为

k_eff = k + (γ/(βΔt))·c + (1/(βΔt²))·m

每步需解一次方程（单自由度时只是一次除法），计算量略大，但当 γ ≥ 1/2、β ≥ γ/2 时方法无条件稳定，步长只受精度约束而不受稳定性约束。工程中最常用的是平均加速度法（β = 1/4，γ = 1/2），它二阶精度、无人工阻尼；若希望抑制高频数值振荡，可改用 HHT-α 或取 γ 略大于 1/2 引入数值阻尼。

稳定性与精度的权衡

稳定性回答“解会不会发散”，精度回答“解离真值有多远”。显式方法受 Δt_cr 硬约束，但每步便宜；隐式方法步长自由，但每步要解方程，且过大的步长会引入周期延长与幅值衰减的数值误差。对单自由度系统，一个实用的判断流程是：先算 ω = √(k/m)，若采用显式，取 Δt ≤ 0.1·(2/ω) 留出安全裕度；若采用隐式，取 Δt ≤ T/20（T = 2π/ω 为固有周期）以保证波形分辨率。本例 T ≈ 0.628 s，隐式取 Δt = 0.03 s 左右即可较好还原波形。

精度还体现在能量行为上。无阻尼自由振动时，显式中心差分的数值能量在稳定域内近似守恒但略有漂移；隐式平均加速度法在无阻尼时同样保持能量，但若取 γ > 1/2，则会人为耗散能量，表现为幅值缓慢衰减——这不是物理阻尼，判读结果时必须区分。

案例一：自由振动

取 m = 1 kg、k = 100 N/m、c = 0.2 N·s/m（ξ ≈ 0.01），初始位移 u₀ = 0.05 m、初速度为零。解析解为 u(t) = 0.05·e^−0.1t·cos(ω_dt)，ω_d ≈ 9.9995 rad/s。用显式中心差分取 Δt = 0.02 s 推进 5 s，可观察到位移峰值从 0.05 m 缓慢衰减到约 0.030 m，与解析包络 e^−0.1t 吻合；若把 Δt 改为 0.25 s，位移峰值会逐周期增大，说明已越过稳定边界。用隐式平均加速度法取 Δt = 0.05 s，波形与解析解几乎重合；若改用 γ = 0.6 的算法阻尼版本，则会看到额外衰减，此时应意识到这是数值耗散而非物理阻尼。

案例二：受迫振动与共振

在同一系统上施加简谐激励 F(t) = F₀·sin(Ωt)，取 F₀ = 1 N、Ω = 10 rad/s，恰好接近固有频率，进入共振区。稳态幅值约为 F₀/(c·ω) = 1/(0.2×10) = 0.5 m，远大于静位移 F₀/k = 0.01 m，放大系数约 50。数值模拟时，显式方法需保证 Δt 远小于激励周期 2π/Ω ≈ 0.628 s，取 Δt = 0.01 s 可分辨每个周期约 60 个点；隐式方法取 Δt = 0.02 s 也能稳定推进。判读要点是：先让瞬态衰减（约 5/ξ 个周期，本例约 500 个周期，实际可只观察前若干周期趋势），再看稳态幅值是否收敛到理论值附近；若幅值持续漂移或出现高频毛刺，应检查步长与算法阻尼设置。

把上述流程固化为可复用的检查顺序：确定 m、c、k 与激励 → 计算 ω、T、ξ → 显式核算 Δt_cr 并留裕度，隐式按 T/20 取步长 → 推进求解并记录位移时程 → 与解析包络或稳态理论值对比 → 若发散查步长，若衰减异常查算法阻尼，若幅值偏低查步长与共振分辨。掌握这一闭环，就具备了对更复杂多自由度与连续体模型进行时间积分选型的基础判断力。