别再死记硬背了！手把手带你推导Newmark-β法（线性加速度法）的核心递推公式-平芜编程栈

从零推导Newmark-β法：结构动力学中的时间积分艺术

想象你正面对一座由微分方程构成的高墙——结构动力学中的运动方程。教科书上跳跃的推导步骤和神秘的积分符号让人望而生畏，而Newmark-β法作为经典的时间积分方法，其核心公式常常被直接给出，背后的推导逻辑却鲜有详细解释。本文将带你一步步拆解这个"黑箱"，用工程师的思维而非数学家的抽象，还原线性加速度假设下递推公式的完整诞生过程。无论你是正在啃教科书的学生，还是需要温故知新的工程师，这里没有"显然可得"的跳跃，只有手把手的地毯式推导。

1. 基础准备：理解动力方程的离散化需求

结构动力学的基本运动方程可以表示为：

m\ddot{u} + c\dot{u} + ku = f(t)

这个看似简洁的方程却包含了位移u对时间t的二阶导数，解析求解仅在极简单情况下可行。数值分析中，我们需要将连续的时间离散化为若干时间步，在每个时间步内对加速度、速度和位移的变化做出合理假设，这就是时间积分法的核心思想。

为什么选择线性加速度假设？在众多可能的假设中，线性加速度（即假定加速度在时间步长Δt内呈线性变化）提供了良好的平衡：

比恒定加速度假设更精确
比高阶多项式假设更简单稳定
物理意义直观，易于实现

下表对比了几种常见的时间积分假设：

假设类型	加速度变化	计算复杂度	精度阶次
恒定加速度	不随时间变化	低	O(Δt²)
线性加速度	随时间线性变化	中	O(Δt³)
三次多项式	随时间二次变化	高	O(Δt⁴)
精确积分	遵循精确动力学方程	极高	精确

提示：Newmark-β法中，当参数β=1/6时，对应线性加速度法，这也是本文推导的特定情况。

2. 线性加速度假设的数学表达

我们明确线性加速度假设的数学表述：在时间区间[ti, ti+1]内，加速度随时间线性变化，即：

\ddot{u}(t) = \ddot{u}_i + \frac{\ddot{u}_{i+1} - \ddot{u}_i}{\Delta t}(t - t_i)

这个看似简单的假设是整个推导的基石。为了从加速度得到速度和位移，需要进行两次积分运算，这正是大多数教材容易跳跃的关键步骤。

积分前的准备工作：

确认积分区间：从ti到t，其中t∈[ti, ti+1]
识别常数项和变量项：式中ui和u˙i是常数，而(t-ti)是变量
准备积分基本公式：
- ∫(t-ti)dt = (1/2)(t-ti)²
- ∫(t-ti)²dt = (1/3)(t-ti)³

3. 第一次积分：从加速度到速度

根据基本微积分，速度是加速度的积分。对线性加速度表达式进行第一次积分：

\dot{u}(t) = \int_{t_i}^t \ddot{u}(\tau) d\tau + C_1

这里有几个关键细节需要注意：

积分变量从t改为τ以避免混淆
必须添加积分常数C1
初始条件决定常数：当t=ti时，u˙(t)=u˙i

展开积分运算：

\begin{aligned} \dot{u}(t) &= \int_{t_i}^t \left[ \ddot{u}_i + \frac{\ddot{u}_{i+1} - \ddot{u}_i}{\Delta t}(\tau - t_i) \right] d\tau + \dot{u}_i \\ &= \ddot{u}_i(t - t_i) + \frac{\ddot{u}_{i+1} - \ddot{u}_i}{2\Delta t}(t - t_i)^2 + \dot{u}_i \end{aligned}

注意：红色标注的u˙i是积分常数，在实际计算中容易被遗漏。它不是被积函数的一部分，而是由初始条件确定的固定值。

4. 第二次积分：从速度到位移

同样地，位移是速度的积分。对得到的速度表达式进行第二次积分：

u(t) = \int_{t_i}^t \dot{u}(\tau) d\tau + C_2

同样需要考虑：

新的积分常数C2
初始条件：当t=ti时，u(t)=ui

展开计算过程：

\begin{aligned} u(t) &= \int_{t_i}^t \left[ \dot{u}_i + \ddot{u}_i(\tau - t_i) + \frac{\ddot{u}_{i+1} - \ddot{u}_i}{2\Delta t}(\tau - t_i)^2 \right] d\tau + u_i \\ &= u_i + \dot{u}_i(t - t_i) + \frac{1}{2}\ddot{u}_i(t - t_i)^2 + \frac{\ddot{u}_{i+1} - \ddot{u}_i}{6\Delta t}(t - t_i)^3 \end{aligned}

常见错误警示：

遗漏积分常数ui（红色标注项）
混淆u¨i和(u¨i+1 - u¨i)/Δt的角色：前者是系数，后者是斜率
积分时错误地将所有项视为变量，实际上u˙i和u¨i在当前步是已知常数

5. 时间步递推公式的最终形式

现在，我们将t取为ti+1（即t = ti + Δt），得到递推公式：

速度更新公式：

\dot{u}_{i+1} = \dot{u}_i + \frac{\ddot{u}_i + \ddot{u}_{i+1}}{2} \Delta t

位移更新公式：

u_{i+1} = u_i + \dot{u}_i \Delta t + \frac{2\ddot{u}_i + \ddot{u}_{i+1}}{6} \Delta t^2

为便于编程实现，我们通常将其重写为更紧凑的形式：

# Python伪代码实现Newmark-β法(β=1/6)的核心更新 def newmark_step(u_i, u_dot_i, u_ddot_i, u_ddot_ip1, dt): u_dot_ip1 = u_dot_i + 0.5 * (u_ddot_i + u_ddot_ip1) * dt u_ip1 = u_i + u_dot_i * dt + (2*u_ddot_i + u_ddot_ip1) * dt**2 /6 return u_ip1, u_dot_ip1

6. 验证与数值稳定性讨论

为了验证我们的推导，考虑一个简单单自由度系统：m=1, c=0, k=4，初始条件u(0)=1, u˙(0)=0，外力f(t)=0。解析解应为u(t)=cos(2t)。

使用Δt=0.1s的线性加速度法数值解与解析解对比如下：

时间(s)	数值解	解析解
0.1	0.98007	0.98007
0.2	0.92106	0.92106
0.5	0.54030	0.54030
1.0	-0.41615	-0.41615

注意：此特例中线性加速度法给出了精确解，这是因为问题本身和解的特性导致的特殊情形。一般情况下，数值解会有一定误差。

关于数值稳定性，线性加速度法（β=1/6）的条件稳定，其稳定条件为：

\Delta t \leq \frac{T_{\text{min}}}{\pi}

其中Tmin是系统的最小周期。在实际应用中，建议取Δt ≤ Tmin/10以保证精度。

7. 工程应用中的实用技巧

在实际结构分析中，线性加速度法的实现还需要考虑以下实际问题：

初始加速度计算：

\ddot{u}_0 = \frac{f(0) - c\dot{u}_0 - ku_0}{m}

迭代求解流程：
- 已知：ui, u˙i, u¨i
- 求解：u¨i+1 通过运动方程
- 更新：ui+1, u˙i+1
- 重复

变步长策略：

# 自适应步长控制示例 def adaptive_step(u, u_dot, u_ddot, dt, tolerance): while True: u_new, u_dot_new = newmark_step(u, u_dot, u_ddot, dt) error = estimate_error(u, u_new) if error < tolerance: return u_new, u_dot_new else: dt *= 0.8 # 减小步长