📚 课程目录第2章 LS-DYNA爆炸分析基础操作

2. LS-DYNA爆炸分析基础操作

爆炸载荷下的结构响应分析，本质上是求解一个强非线性、大变形、含接触与材料失效的瞬态动力学问题。LS-DYNA 采用显式中心差分法推进时间步，理解它的求解逻辑，是正确设置关键字、判断结果是否可信的前提。本节围绕显式动力学求解逻辑、单位制选择、关键字文件结构、前处理与后处理工具链，以及一个可跑通的最小算例展开。

一、显式动力学求解逻辑

显式方法求解的运动方程为 M·a = F_ext − F_int，其中 M 为集中质量矩阵。由于采用集中质量（对角质量矩阵），无需组装和求逆整体刚度矩阵，加速度可直接由节点力除以节点质量得到：a = (F_ext − F_int)/M。随后用中心差分更新速度与位移：

v(t+Δt/2) = v(t−Δt/2) + a(t)·Δt
d(t+Δt) = d(t) + v(t+Δt/2)·Δt

这套递推是条件稳定的，稳定条件由 Courant 准则给出：Δt ≤ Δt_crit = 2/ω_max，工程上常写成 Δt = α·L_min/c，其中 L_min 是最小单元特征长度，c 是材料声速，α 为时间步缩放因子（默认约 0.9）。这意味着：网格中最小的一两个单元会决定整个模型的求解步长。若某处存在 0.5 mm 的碎小单元，而主体网格为 5 mm，则时间步会被拉小一个数量级，计算时间成倍增加。因此爆炸分析中要严格控制最小单元尺寸，必要时用 *CONTROL_TIMESTEP 的质量缩放（DT2MS）来人为放大步长，但必须监控质量增加比例（一般控制在 1%~5% 以内）。

显式方法的另一特点是"一步一算、不迭代"，单步成本低但步数极多。爆炸载荷持续时间通常在毫秒量级，而稳定步长常在微秒甚至纳秒量级，因此总步数可达十万到百万级。这也解释了为什么爆炸分析中要尽量使用单点积分实体单元、避免过多小单元、并合理设置终止时间。

二、单位制选择

LS-DYNA 本身不带单位，所有输入必须自洽。爆炸分析中最常用的是 mm–ms–kg–GPa 与 m–s–kg–Pa 两套。推荐使用 mm–ms–kg–GPa：长度 mm、时间 ms、质量 kg、力 kN、应力 GPa、密度 kg/mm³（如钢 7.85e-6 kg/mm³）、能量 kJ。这套单位下，钢的声速约 5.1 mm/ms，与 mm 网格配合时步长数值直观。

常见错误是把密度写成 7850 kg/m³ 却用 mm 建模，导致质量放大 10⁹ 倍，结构几乎不动。检查方法：用 *DATABASE_ASCII 输出总质量，与手算体积×密度对比。另一处易错是重力加速度，在 mm–ms 制下 g = 9.81e-3 mm/ms²，若忽略重力对爆炸这种毫秒级过程影响不大，但涉及准静态阶段时不可漏。

三、关键字文件结构

一个完整的 K 文件由若干关键字卡片块组成，顺序不严格但逻辑上分为：控制卡、材料卡、部件与单元、接触、载荷、输出。最小爆炸算例通常包含：

*KEYWORD
*CONTROL_TERMINATION
$  ENDTIM 终止时间(ms)
   0.5
*CONTROL_TIMESTEP
$  DTINIT  TSSFAC  DT2MS
   0.0     0.9     0.0
*CONTROL_ENERGY
$  HGEN  RWEN  SLNTEN  RYLEN
   2     2     2      2
*DATABASE_BINARY_D3PLOT
$  DT 输出间隔(ms)
   0.01
*PART
$  PID  SECID  MID  EOSID  HGID
   1    1      1    1      0
*SECTION_SOLID
$  SECID  ELFORM
   1      1
*MAT_PLASTIC_KINEMATIC
$  MID  RO(kg/mm3)  E(GPa)  PR  SIGY(GPa)  ETAN  BETA
   1    7.85e-6     210     0.3 0.3         2.0   0.0
*EOS_JWL
$  EOSID  A  B  R1  R2  OMEG  E0  V0
   1      371 3.23 4.15 0.95 0.30 7.0 1.0
*INITIAL_VOLUME_FRACTION_GEOMETRY
...
*CONTACT_AUTOMATIC_SURFACE_TO_SURFACE
...
*SET_NODE_LIST
...
*LOAD_BLAST_ENHANCED
$  BLAST  ...  TNT当量(kg)  起爆点坐标
...
*END

注意 *MAT_PLASTIC_KINEMATIC 中 ETAN 为切线模量，BETA 为各向同性硬化参数（0 为随动硬化，1 为各向同性硬化）。爆炸下结构常经历反向加载，BETA 取 0 更接近实际。JWL 参数需与炸药类型匹配，TNT 常用 A=371 GPa、B=3.23 GPa、R1=4.15、R2=0.95、ω=0.30、E0=7.0 GPa。

四、前处理与后处理工具链

前处理常用 LS-PrePost（免费、原生支持 K 文件）、HyperMesh、ANSYS Workbench 的 Explicit Dynamics 模块。LS-PrePost 可直接读入 K 文件、检查关键字拼写、显示单元质量与最小尺寸。后处理同样以 LS-PrePost 为主，读取 d3plot 查看应力云图、位移时程、能量曲线。判断结果是否可信，重点看三件事：沙漏能（Hourglass Energy）应小于总内能的 5%；总能量应基本守恒；质量增加应控制在 5% 以内。若沙漏能异常大，说明单点积分单元出现零能模式，需加密网格或改用 ELFORM=2（全积分）或添加 *CONTROL_HOURGLASS。

五、最小算例跑通

以一个 100 mm × 100 mm × 5 mm 的钢板为例，四边固支，中心正上方 200 mm 处起爆 0.5 kg TNT。网格 2 mm 六面体，单元数约 6250。设置终止时间 0.5 ms，D3PLOT 每 0.01 ms 输出一帧。提交计算后，先看 d3plot 第一帧确认模型与边界正确，再看板中心位移时程是否在 0.1~0.2 ms 达到峰值（示例中约 8~12 mm 量级），最后检查能量曲线。若计算在几秒内结束且位移为零，多半是单位制或载荷卡片错误；若时间步极小、计算停滞，检查最小单元尺寸与 DT2MS 设置。跑通这个算例后，再逐步加入应变率效应（*MAT_PLASTIC_KINEMATIC 的 C、P 参数）、失效准则与接触，即可过渡到真实工程问题。