📚 课程目录第2章显式动力学基础与LS-DYNA求解原理

2. 显式动力学基础与LS-DYNA求解原理

在车辆碰撞安全仿真中，物理过程通常只持续 80～150 ms，却伴随大变形、接触、材料失效与断裂。若采用隐式方法求解，每一时间步都要组装并求解大型刚度矩阵，还要反复迭代判断收敛，代价极高且容易发散。LS-DYNA 的核心是显式时间积分：它不求解联立方程组，而是用当前时刻已知的量直接推算下一时刻的量，因此每个时间步计算量小、天然适合并行，也能稳定穿过接触状态突变和单元畸变。

中心差分：显式积分的基本形式

LS-DYNA 默认采用中心差分法（central difference）。以节点位移 u、速度 v、加速度 a 为例，其递推关系为：

a(t) = M⁻¹ · (F_ext(t) − F_int(t))
v(t+Δt/2) = v(t−Δt/2) + a(t)·Δt
u(t+Δt) = u(t) + v(t+Δt/2)·Δt

其中 M 为质量矩阵，F_ext 为外力（碰撞中主要是接触力与惯性力），F_int 为单元内力。关键点在于：LS-DYNA 使用集中质量矩阵（lumped mass），M 是对角阵，求逆只是逐节点做除法，无需矩阵分解。这就是显式方法“每步便宜”的根本原因。速度定义在半个时间步上（交错网格），位移和加速度定义在整步上，这种交错存储保证了二阶精度。

需要特别注意：显式方法有条件稳定，即时间步长不能任意取大。若 Δt 超过临界值，数值解会指数发散，表现为能量异常增长、单元飞散。因此时间步长控制是显式求解的第一生命线。

条件稳定与临界时间步长

临界时间步长由结构中最小单元与材料波速共同决定。对杆单元，其形式为 Δt_crit = L/c，其中 L 为单元特征长度，c 为材料声速。对三维实体单元，LS-DYNA 采用更保守的估算：

Δt_crit = L_s / c
c = sqrt(E / ρ)          （弹性波速，E 为弹性模量，ρ 为密度）
L_s = V / A_max          （特征长度，V 为单元体积，A_max 为最大面面积）

实际计算中，LS-DYNA 会扫描全部单元取最小值，再乘以安全系数（默认 0.9，可用 *CONTROL_TIMESTEP 的 TSSFAC 修改）。举例：某钢制薄壁吸能盒单元特征长度 L_s = 3 mm，E = 210 GPa，ρ = 7.85×10⁻⁹ t/mm³，则 c ≈ 5170 mm/ms，Δt_crit ≈ 5.8×10⁻⁴ ms。若整车模型最小单元为 1 mm，Δt 会降到约 1.9×10⁻⁴ ms，求解 100 ms 需要约 52 万步，计算量急剧上升。

单位制必须自洽。车辆碰撞常用 mm–ms–t–N–MPa 体系：长度 mm、时间 ms、质量 t、力 N、应力 MPa、密度 t/mm³。钢的密度写作 7.85e-9，而不是 7850。

质量缩放：用可控误差换计算效率

当模型中存在少量极小单元时，它们会把全局时间步长拖到很低。质量缩放（mass scaling）的思路是：人为增大这些单元的质量密度，使波速 c 下降、临界步长上升。由 Δt ∝ sqrt(ρ)，若希望把 Δt 提高 k 倍，密度需增大 k² 倍。

LS-DYNA 提供两类质量缩放：

固定时间步质量缩放：在 *CONTROL_TIMESTEP 中设置 DT2MS 为目标步长（负值表示只对步长小于该值的单元缩放，且不改变总质量分布逻辑）。例如 DT2MS = −1.0e-3 ms，则所有 Δt 小于 1.0e-3 ms 的单元被自动加质量。
选择性质量缩放：通过 *PART 或 *ELEMENT_MASS 对特定部件、特定单元单独处理，避免污染关键结构。

质量缩放不是免费的。增加的质量会改变惯性效应，若总质量增加超过 1%～2%，碰撞响应（尤其是加速度峰值和吸能分配）可能明显失真。因此必须监控 glstat 中的 “added mass” 与 “total mass” 比例，并检查 matsum 中动能是否异常。工程经验是：整车碰撞模型质量增加控制在 1% 以内，部件级模型可放宽到 3%～5%，但需做无缩放对比验证。

沙漏控制：单点积分的必要代价

为追求速度，LS-DYNA 实体单元默认采用单点积分（reduced integration）。单点积分只有一个积分点，计算量最小，但会产生零能模式——即沙漏（hourglass）。沙漏表现为单元在无外力做功的情况下发生交替畸变，网格呈“棋盘”状，能量凭空产生，最终导致结果完全不可信。

沙漏控制的本质是给单元施加一个抵抗非物理变形模式的刚度或黏性力。常用方法包括：

黏性沙漏控制（type 1、2、3）：默认类型，对小变形有效，但可能引入过多阻尼。
刚度沙漏控制（type 4、5）：适用于大变形和金属塑性问题，对碰撞仿真更稳健。
全积分单元：如 ELFORM=2 的实体单元，无沙漏但计算量成倍增加，且可能发生剪切锁定。

在 *HOURGLASS 关键字中可指定类型与系数 QH。对薄壁吸能结构，推荐刚度型沙漏控制，QH 取 0.03～0.10。判断沙漏是否受控，应查看 glstat 中 “hourglass energy” 与总内能之比：一般要求小于 5%，关键吸能部件应小于 2%。若沙漏能持续增长，说明网格过粗或加载方式激发了零能模式，需要细化网格或调整控制参数。

三个控制量的协同关系

时间步长、质量缩放与沙漏控制并非彼此独立。网格细化会降低 Δt，迫使更多质量缩放；质量缩放过大又改变惯性效应，可能掩盖沙漏问题；沙漏控制过强则引入额外刚度，使结构偏硬、吸能量偏高。因此实际建模中应遵循以下判断顺序：

先检查网格质量，消除长宽比过大、翘曲严重的单元，从源头避免极小 Δt。
再设定 DT2MS，使全局步长提升到可接受范围，同时监控 added mass 比例。
最后根据部件变形模式选择沙漏控制类型与系数，用 glstat 和 matsum 验证能量守恒。

不要用质量缩放去“救”一个网格质量很差的模型。若 added mass 超过 5%，即使计算能跑完，碰撞力曲线和吸能分配也已不可信。正确做法是回到网格划分阶段，控制最小单元尺寸与单元质量指标。

理解上述原理后，后续的接触设置、材料模型与求解控制才有稳定的数值基础。显式求解的每一步都建立在“当前已知量直接推算下一步”之上，而条件稳定、质量缩放与沙漏控制正是保证这条递推链不发散、不失真的三道闸门。