Newmark-β法瞬态结构分析

简单说：Newmark-β法是隐式积分，无条件稳定（时间步可以比较大），适合低频主导的结构振动、地震响应等中长时间历程问题；中心差分法是显式积分，有条件稳定（时间步必须非常小，受最小网格尺寸约束），适合冲击、爆炸等高频短时问题。以汽车碰撞仿真为例：座椅静态强度分析用隐式（Newmark），正面碰撞100ms时程用显式（LS-DYNA中心差分）。

稳定是指计算不会数值发散（误差不会指数级放大），但绝对不等于精度没问题。时间步太大会导致两个问题：①高频成分被数值阻尼过滤掉（人工阻尼）；②低频响应有相位延迟误差。实务中时间步通常取最高关注频率周期的1/10～1/20，在稳定性和计算精度之间取平衡。

Newmark（1959年）假设时间步 $[t_n, t_{n+1}]$ 内加速度按参数 $\beta$ 插值：

积分两次得到位移和速度更新公式：

参数 $\beta$ 控制加速度插值权重（决定稳定性），$\gamma$ 控制速度插值（决定精度和数值阻尼）：

对单自由度系统（SDOF）分析Newmark法的稳定性，将递推格式写成放大矩阵 $\mathbf{A}$：

无条件稳定要求放大矩阵谱半径 $\rho(\mathbf{A}) \le 1$（对所有 $\Omega = \omega\Delta t$）。

对Newmark法，无条件稳定的充要条件为：

平均加速度法（$\beta=1/4, \gamma=1/2$）恰好满足等号，是无条件稳定的最低参数配置。$\gamma > 1/2$ 会引入数值阻尼，降低精度到1阶——这就是为什么HHT-α法从运动方程修改而非改动 $\gamma$ 来引入阻尼。

HHT-α法（Hilber-Hughes-Taylor，1977年）的巧妙之处是在运动方程本身引入加权，而不是改变积分格式：

参数 $-1/3 \le \alpha \le 0$（通常取 $-0.3 \le \alpha \le -0.05$），对应参数：

这样既保持二阶精度（$\gamma=1/2$ 时2阶），又对高频内容提供可控数值阻尼（$|\alpha|$ 越大阻尼越强）。$\alpha = 0$ 时退化为Newmark平均加速度法。HHT-α是目前所有主流FEM软件默认隐式积分算法的基础。

每个时间步 $[t_n, t_{n+1}]$ 的计算流程：

1. 预测（Predictor）： 用已知 $\{u\}_n, \{\dot{u}\}_n, \{\ddot{u}\}_n$ 预测：
   $\{u\}_{n+1}^{(0)} = \{u\}_n + \Delta t\{\dot{u}\}_n + \Delta t^2(\frac{1}{2}-\beta)\{\ddot{u}\}_n$
   $\{\dot{u}\}_{n+1}^{(0)} = \{\dot{u}\}_n + (1-\gamma)\Delta t\{\ddot{u}\}_n$
2. 组装等效刚度矩阵： $\hat{K} = K + \frac{\gamma}{\beta\Delta t}C + \frac{1}{\beta\Delta t^2}M$
3. Newton-Raphson迭代（非线性时）： 求解 $\hat{K}\Delta u = \hat{f}_{residual}$ 直到收敛
4. 校正（Corrector）： 用 $\Delta u$ 更新位移、速度、加速度
5. 输出应力/应变，推进到下一步

步骤2中 $\hat{K}$ 的组装和LU分解是最耗时的部分。对于线性问题只需分解一次；非线性问题每步至少分解一次，用BFGS近似可减少重分解次数。

非线性问题的几个关键点：

• 自适应时间步： 根据每步Newton-Raphson迭代次数自动调整 $\Delta t$（迭代次数多→减小，顺利收敛→放大）。Abaqus的自动增量控制（AINC）是业界标杆。
• 接触状态突变： 接触开/闭状态改变时需要事件驱动的时间步细化（Event-driven Bisection），否则出现穿透或接触力振荡（spurious oscillation）。轻微数值阻尼（HHT-α, α≈-0.05）能有效抑制接触振荡。
• 刚度更新频率： Full Newton（每迭代重组刚度矩阵）精度最高；Modified Newton（多步用同一切线刚度）节省计算；Quasi-Newton（BFGS）是折中。
• 能量守恒检查： 每步验证 $\Delta E_{kin} + \Delta E_{int} - W_{ext} = 0$（能量残差/输入能量 < 1%），异常增加提示积分发散。

直接积分法（Newmark）不能对每阶模态单独设置阻尼比，只能用Rayleigh阻尼：

对角频率 $\omega_i$ 的模态，等效阻尼比为：

标定方法：选两个目标频率 $\omega_1$（低频，如一阶固有频率）和 $\omega_2$（高频，关注上限），指定两者的阻尼比 $\xi_1 = \xi_2 = \xi$（工程结构通常2～5%），解方程组：

注意：频率范围以外的模态阻尼比会偏高（特别是高频），$\beta_R K$ 项对高频模态引入很大阻尼——这就是为什么Rayleigh阻尼会"吸收"高频噪声。

地震响应分析的典型流程：

1. 模态分析预置： 先确认结构前几阶固有频率（对应周期是否在地震影响范围0.1～3秒内）
2. 时间步确定： $\Delta t \le T_{min}/20$。关注33 Hz（混凝土结构建筑高频限）时 $\Delta t \le 0.03/20 = 0.0015\,\text{s}$；地震波本身通常以0.005 s或0.01 s采样，计算时间步应不大于采样间隔。
3. 阻尼设置： 钢结构典型 $\xi = 2\%$，混凝土 $\xi = 5\%$。Rayleigh阻尼标定频率取第1阶和主要响应模态。
4. 输入地震波： 实测加速度时程（如Kobe 1995）或规范人工波，施加为基础加速度荷载
5. 结果提取： 楼层最大加速度、层间位移角（Interstory Drift Ratio）、构件应力时程

高频振荡是平均加速度法（$\alpha=0$）无数值阻尼的表现。解决方案：

• 启用HHT-α数值阻尼（$\alpha = -0.05$～$-0.1$），抑制高于关注频率的分量
• 检查时间步是否足够小（高频成分被欠采样会产生混叠振荡）
• 后处理时用低通滤波（截止频率取关注最高频率的2倍）
• 若是接触冲击引起的振荡，检查接触刚度（惩罚法的惩罚系数过大会引起高频振荡）

主要软件的瞬态积分方案对比：

Abaqus用 $\rho_\infty$（无穷频率处的谱半径，0≤$\rho_\infty$≤1）参数化：

对应关系：$\rho_\infty = 1$ → $\alpha = 0$（无阻尼，Newmark平均加速度）；$\rho_\infty = 0$ → $\alpha = -1/3$（最大阻尼，HHT-α极端情况）。Abaqus默认 $\rho_\infty = 0.9929$（对应 $\alpha \approx -0.005$，极轻数值阻尼）。一般非线性接触问题推荐 $\rho_\infty = 0.9$ 左右（$\alpha \approx -0.05$）。

几个值得关注的方向：

• GSSSS统一框架（Generalized Single Step Single Solve）： 将Newmark、HHT-α、WBZ等算法统一在一个广义框架下，参数化空间更大，可在精度、稳定性和数值阻尼之间更自由地权衡，实现用户指定的频率响应函数。
• 能量动量守恒算法（Energy-Momentum Method）： 在大位移几何非线性问题（柔性体、生物力学）中无条件稳定且严格守恒，消除了Newmark在大位移时的能量漂移问题。
• 基于机器学习的自适应时间步： 用LSTM预测下一步非线性程度，提前调整 $\Delta t$，减少不必要的Newton迭代。
• IGA（等几何分析）+ 时间积分： B样条/NURBS基函数在时间维度的应用，实现空间和时间都达到高阶（$p$阶）精度，对光滑动力学问题误差远小于传统FEM+Newmark组合。
• 显隐式耦合（Operator Splitting）： 对同一模型中的"硬"区域（细网格、高刚度）用隐式，"软"区域用显式，各自用最优时间步，节省计算量。

两种方法的核心对比：

选择原则：线性结构+频率响应分析→模态叠加（Nastran SOL 111/112）；含接触、材料非线性、大变形→直接积分（Abaqus *DYNAMIC, Ansys FULL法）。