📚 课程目录第2章显式动力学理论基础

2. 显式动力学理论基础

跌落冲击仿真的本质，是求解结构在极短时间（通常 0.5～20 ms）内承受大位移、大应变、接触碰撞的动力学响应。Ansys Mechanical 中的显式动力学求解器（如 LS-DYNA、Autodyn 或 Explicit Dynamics 模块）与静力学、隐式动力学在数学结构上有根本区别。理解运动方程、中心差分法、时间步长、沙漏控制与质量缩放这五个概念，是判断一次跌落仿真结果是否可信的基础。

一、运动方程：从虚功原理到半离散形式

显式动力学求解的是动力学平衡方程。对连续体，用虚功原理或加权残量法离散后，得到半离散运动方程：

M·a(t) + C·v(t) + F_int(u,t) = F_ext(t)

其中：
M      质量矩阵（集中质量，对角阵）
C      阻尼矩阵
F_int  内力向量（单元应力经形函数积分得到）
F_ext  外力向量（重力、接触力、压力等）
a, v, u 节点加速度、速度、位移

与隐式方法不同，显式方法采用集中质量矩阵（lumped mass），即把单元质量按节点分配，使 M 成为对角阵。这样 M 求逆只需逐元素取倒数，无需矩阵分解，这是显式方法每步计算量极小的关键。跌落仿真中，重力加速度 g=9806.65 mm/s²（或 9.80665 m/s²）作为体积力施加，接触力通过罚函数或运动约束算法在 F_ext 中体现。

二、中心差分法：时间推进的核心

中心差分法（Central Difference Method）是显式求解的时间积分格式。它用相邻两个时刻的位移差表示速度，用相邻两个时刻的速度差表示加速度：

v(t) = [u(t+Δt) - u(t-Δt)] / (2Δt)
a(t) = [u(t+Δt) - 2u(t) + u(t-Δt)] / Δt²

代入运动方程并整理，得到位移递推公式：

u(t+Δt) = 2u(t) - u(t-Δt) + Δt²·M⁻¹·[F_ext(t) - F_int(t) - C·v(t)]

该式表明：下一时刻位移只依赖当前和上一时刻的已知量，无需迭代求解方程组，因此称为“显式”。每一步计算量正比于单元数，适合百万级单元的大模型。但中心差分法是有条件稳定的：时间步长必须小于临界值，否则位移会指数发散。

三、时间步长与临界时间步

临界时间步由系统中最小单元的特征尺寸和材料波速决定。对一维杆单元，稳定条件为：

Δt_crit = L_e / c
c = sqrt(E/ρ)   （弹性波速）

其中：
L_e  单元特征长度（最小边长或面积/最长边）
E    弹性模量
ρ    密度

对三维实体单元，常用保守估计：

Δt_crit = L_e / [c·(1 + 泊松比影响系数)] ≈ 0.8～0.9 × L_e/c

实际求解器会扫描所有单元，取全局最小值作为时间步，并乘以安全系数（默认 0.9）。例如：钢材 E=210 GPa，ρ=7850 kg/m³，c≈5170 m/s；若最小单元尺寸 1 mm，则 Δt_crit≈0.19 μs。要模拟 5 ms 跌落，需要约 26000 步。若网格中有一个 0.1 mm 的微小单元，步长会骤降到 0.019 μs，步数增加 10 倍——这就是“一个坏单元拖慢整个模型”的原因。

四、沙漏控制：单点积分的代价

显式求解为降低计算量，常采用单点积分（reduced integration）实体单元。单点积分只在一个积分点计算应力，效率高，但会产生“沙漏模式”——单元在无应变能增加的情况下发生交替变形，表现为网格呈锯齿状、棋盘状振荡。这种零能模式若不控制，会污染位移场和应力结果。

Ansys/LS-DYNA 提供两类沙漏控制：

刚度控制（IHQ=4 等）：对沙漏模式施加弹性刚度，适合冲击初期；
粘性控制（IHQ=1～3）：施加与变形速率相关的阻尼，适合高速碰撞。

工程判断标准：沙漏能应小于总能量的 5%。若超过，说明网格过粗或控制系数过小，应优先细化网格，而不是盲目增大沙漏系数（过大系数会使响应偏刚）。跌落仿真中，塑料外壳、薄壁件最容易出现沙漏，建议对这些区域使用全积分单元或至少 2 层以上网格。

五、质量缩放：用精度换时间

质量缩放（Mass Scaling）通过人为增大密度 ρ 来提高波速分母中的质量项，从而增大临界时间步。由 Δt_crit ∝ √(ρ/E)·L_e 可知，密度增大 100 倍，步长增大 10 倍，步数减少 10 倍。但质量增加会改变惯性效应，使冲击响应偏慢、峰值力偏低。

质量缩放分两种：

全局缩放：对所有单元统一增大密度，简单但可能严重失真；
局部/选择性缩放：只对小于目标步长的单元缩放，控制质量增加总量。

在 LS-DYNA 中可用 *CONTROL_TIMESTEP 设置 DT2MS（目标时间步），负值表示只对小于该步长的单元缩放。例如设置 DT2MS = −1e−6 s，则所有 Δt_crit < 1 μs 的单元被缩放至 1 μs。求解结束后必须检查 glstat 文件中的 “mass increase” 和 “added mass” 百分比。

跌落仿真经验值：质量增加控制在 1% 以内对冲击峰值影响通常可接受；若达到 5%～10%，需对比无缩放模型验证峰值加速度和接触力偏差。对于手机、路由器等消费电子跌落，建议先做无缩放基准，再逐步放宽 DT2MS。

六、五个概念的联动关系

这五个概念不是孤立的：运动方程决定求解对象；中心差分法决定时间推进方式；临界时间步决定步长上限；沙漏控制保证单点积分单元不产生零能振荡；质量缩放则在步长过小时提供“合法加速”手段。一次可信的跌落仿真，应当同时满足：

时间步长 ≤ 0.9 × 全局最小 Δt_crit；
沙漏能 < 总能量 5%；
质量增加 < 1%（冲击主导）或 < 5%（含准静态段）；
动能与内能曲线无异常尖峰或负值；
接触力、加速度峰值与网格密度收敛性检查通过。

实际建模时，建议先以较粗网格和较小质量缩放快速试算，观察能量曲线和沙漏能占比；确认稳定后，再细化关键区域网格并收紧质量缩放。这样既避免“一上来就百万单元跑不动”，也避免“为了跑得快而把物理做丢”。