落体冲击试验仿真
理论与物理
自由落体力学
老师，我想用FEM预测手机掉落时的冲击，能先教我一些关于掉落的基本物理知识吗？
问得好。跌落测试的起点是忽略空气阻力的自由落体。从高度 $h$ 处轻轻释放的物体撞击地面瞬间的速度，可以从能量守恒定律推导出来：
例如，将智能手机从 1.2 m 的高度（大约从口袋取出时的高度）掉落：
将这个速度作为“初速度”赋予产品，是跌落FEM的基本设置方法。质量 $m$ 虽然像约定俗成一样消掉了，但当然会影响碰撞后的变形和应力。
4.85 m/s，换算成时速大约是17 km/h吧。想象一下以自行车般的速度撞向地面，感觉还挺快的…。
是的，而且接触时间极短。智能手机的角撞击混凝土地面时，接触时间仅有 0.3〜0.5 ms。动量的变化 $\Delta p = mv$ 在这极短时间内发生，因此平均冲击力为：
对于200 g的物体产生约2,400 N的力，这意味着加速度达到约1,200 G。这就是破坏电子元件焊点连接的原因。
碰撞接触力学
碰撞瞬间，接触力是用什么模型计算的呢？可以简单地像弹簧那样考虑吗？
实际上，对于弹性碰撞，赫兹接触理论是基础。两个弹性体接触时，接触力 $F$ 与接近量（压入量）$\delta$ 的关系为：
这里 $K_H$ 是赫兹接触刚度，由两个物体的弹性模量、泊松比和曲率半径决定。但是，跌落测试中会发生塑性变形或破坏，仅靠赫兹理论是不够的。在FEM中，通过在单元层面求解弹塑性本构关系，可以直接计算非线性接触力。
FEM的接触处理具体是怎么做的呢？
主要有罚函数法和拉格朗日乘子法两种。在实际的跌落分析中，罚函数法占绝大多数。其机制很简单，当两个面试图穿透时，虚拟的弹簧会将其推回：
$k_p$ 是罚刚度，$g_n$ 是间隙（负值表示穿透量）。$k_p$ 太小会导致产品“穿过”地面，太大则会使时间步长变得极小，计算成本爆炸式增长。先用默认值运行，再进行微调是实用的方法。
能量收支与吸收机制
跌落的能量都去哪里了呢？即使没有损坏，有时也不会反弹。
这是非常关键的一点。碰撞前的总能量是动能 $E_k = \frac{1}{2}mv^2$，碰撞后这些能量会分配到多种形态中：
- 弹性应变能 $E_{elastic}$ — 变形后恢复的部分。反弹的能量来源
- 塑性耗散能 $E_{plastic}$ — 消耗于永久变形。外壳凹陷或开裂
- 摩擦能 $E_{friction}$ — 接触面滑动消耗
- 阻尼能 $E_{damping}$ — 材料粘性吸收
- 反弹能 $E_{rebound}$ — 反弹后的动能
用恢复系数（COR）$e$ 表示则为 $E_{rebound} = e^2 \cdot E_k$。智能手机玻璃面撞击瓷砖时，$e \approx 0.2\text{〜}0.4$，大部分能量被塑性变形和内部阻尼消耗。
跌落姿态与应力集中
听说跌落角度不同，结果会大不相同，哪种姿态最严酷呢？
跌落姿态会戏剧性地改变结果。同样是1.2 m跌落，面跌落和角跌落时局部应力可能相差10倍以上：
| 跌落姿态 | 初始接触面积 | 峰值加速度 | 破坏风险 | 典型损伤 |
|---|---|---|---|---|
| 面跌落（平面） | 大（整个面） | 200〜500 G | 低 | 液晶屏破裂、基板弯曲 |
| 边跌落（边缘） | 中（线状） | 500〜1,500 G | 中 | 框架变形、按键破损 |
| 角跌落（角落） | 极小（点状） | 1,000〜5,000 G | 最大 | 外壳开裂、IC剥离 |
IEC 60068-2-31要求进行6个面、12条边、8个角，共计26个方向的跌落。在FEM中，至少也必须进行面、边、角这三种代表性姿态的分析，要考察设计余量的话，理想情况是运行全部26个方向。
26个方向，就算一次分析需要30分钟，也要13小时…这现实吗？
所以高性能计算集群就派上用场了。像三星和苹果这样的大公司，会使用数十台工作站并行运行26个方向，一天内就能得出全部结果。对于中小企业来说，首先聚焦于3个角跌落方向+1个面跌落方向，共计4个案例，以确定最严酷的条件，这是比较现实的。
应变率效应
教科书上写着冲击问题中应变率的影响很大，跌落测试也需要考虑吗？
必须考虑。跌落冲击中典型的应变率是 $10^1 \sim 10^3$ /s，是准静态测试（$10^{-3}$ /s）的一万倍以上。许多金属和塑料在高应变率下屈服应力会上升。代表性的模型是Cowper-Symonds准则：
$\sigma_0$ 是准静态屈服应力，$D$ 和 $q$ 是材料常数。例如聚碳酸酯（PC）的 $D = 10$ /s、$q = 2$ 左右，应变率 100 /s 时屈服应力会跃升至准静态的约3倍。如果忽略这一点，仿真中变形会比实际更大，有误判为“损坏”的风险。
iPhone跌落测试秘闻
据传，苹果在开发初代iPhone时，每天要跌落数十台原型机。iPhone 12之后，由于采用了Ceramic Shield玻璃，据说跌落耐受性提高了4倍，但这“4倍”的背后，是使用LS-DYNA进行的数千个案例的参数化分析。通过用内聚单元模拟玻璃的微裂纹扩展，并以0.1mm为单位优化边缘倒角的曲率，才得到了这个结果。
跌落冲击的控制方程
- 运动方程：$M\ddot{u} + C\dot{u} + F_{int}(u) = F_{ext}(t)$ ── 跌落冲击中惯性项 $M\ddot{u}$ 占主导。加速度可达数千G，因此准静态分析完全不适用
- 接触力：$F_{contact} = k_p \cdot \max(0, -g_n)$ ── 罚函数法产生的接触力。$g_n$ 是间隙距离，仅当其为负（穿透）时才产生力
- 能量守恒：$E_{total} = E_k + E_{internal} + E_{contact} + E_{hourglass}$ ── 总能量随时间变化在初始值的±5%以内是分析健康运行的必要条件
- CFL条件：$\Delta t \leq \frac{L_{min}}{c}$ ── 显式算法的稳定极限。$L_{min}$ 是最小单元长度，$c = \sqrt{E/\rho}$ 是弹性波速。对于钢材 $c \approx 5{,}000$ m/s，1 mm单元则 $\Delta t \leq 0.2$ μs
单位制注意事项（跌落分析中常见）
| 物理量 | SI（米制） | mm-ms-tonne制 | 换算备注 |
|---|---|---|---|
| 长度 | m | mm | 1 m = 1000 mm |
| 时间 | s | ms | 1 s = 1000 ms |
| 质量 | kg | tonne | 1 kg = 0.001 tonne |
| 速度 | m/s | mm/ms | 数值相同（1 m/s = 1 mm/ms） |
| 力 | N | N | 相同 |
| 应力 | Pa | MPa | 1 MPa = 10⁶ Pa |
| 密度 | kg/m³ | tonne/mm³ | 钢: 7.85×10⁻⁹ tonne/mm³ |
| 重力加速度 | 9810 mm/s² | 9.81×10⁻³ mm/ms² | LS-DYNA中使用 mm/ms² |
数值解法与实现
显式算法的时间积分
听说跌落分析几乎总是用显式算法，隐式算法不行吗？
跌落冲击是数毫秒的超短时间现象，接触的发生与消失、大变形、材料塑性化同时发生。隐式算法需要在每个时间步求解联立方程，每当接触条件变化时收敛就会变得困难。基于中心差分法的显式算法则是：
如果将质量矩阵 $M$ 对角化（集中质量），就不需要求解联立方程。可以独立计算每个节点的加速度，因此即使非线性很强的问题也能稳定推进。但这是条件稳定的，稳定的时间步长为：
对于钢材的1 mm六面体单元，$\Delta t \approx 0.2$ μs。要模拟5 ms的碰撞现象，大约需要25,000个时间步。
25,000步，真是惊人的数量。不过每一步计算量轻，所以总体上还是快的…。
没错。智能手机级别的模型（50万个单元），使用LS-DYNA并调用16个核心，30分钟到1小时就能完成。用隐式算法解同样的问题，会因为接触切换导致无法收敛，不知道要花多少天。跌落分析使用显式算法不是因为“快”，而是因为“能解出来”。
接触算法
前辈说接触设置是最难的。具体该怎么做呢？
跌落分析的接触主要分为两大类：
- 产品 vs 地面 — 外部接触。将地面设为刚体面则无需单元，效率高。对于混凝土或瓷砖等非常坚硬的面，刚体近似就足够了
- 产品内部的自接触 — 外壳变形后与内部的基板或电池接触的情况。容易被忽视，但对电子设备很重要
在LS-DYNA中，*CONTACT_AUTOMATIC_SINGLE_SURFACE 是常用选择，它会将所有外表面作为一个接触组自动检测。罚刚度的比例因子（SFS/SFM）通常先用默认值1.0运行，如果发现穿透，再提高到2.0〜5.0，这是实务做法。
LS-DYNA 输入示例
能给我看看具体的LS-DYNA关键字设置吗？假设是从1米高的面跌落。
用mm-ms-tonne制表示，1米跌落的碰撞速度为 $v = \sqrt{2 \times 9.81 \times 1.0} = 4.43$ m/s = 4.43 mm/ms。这样设置：
$ --- 初速度（赋予整个产品集向下的速度）---
*INITIAL_VELOCITY_SET
$ nsid vx vy vz
1, 0.0, 0.0, -4.43
$ --- 重力（仅在追踪反弹时需要）---
*LOAD_BODY_Z
$ lcid sf
1, 1, 9.81e-3 $ mm/ms^2
$ --- 刚体墙（无限坚硬的地面）---
*RIGIDWALL_PLANAR
$ xt yt zt xh yh zh
0., 0., 0., 0., 0., 1.
$ --- 接触（自动单一面）---
*CONTACT_AUTOMATIC_SINGLE_SURFACE
$ ssid msid sstyp mstyp sboxid mboxid spr mpr
1, 0, 3, 0, 0, 0, 1, 1
$ fs fd dc vc
0.3, 0.2, 0., 0.
$ --- 时间控制 ---
*CONTROL_TIMESTEP
$ dtinit tssfac
0., 0.9要点是使用 TSSFAC=0.9，即使用稳定时间步长的90%。默认的0.9通常没问题，但如果接触非常剧烈，有时会降到0.67。
Abaqus/Explicit 输入示例
Abaqus/Explicit怎么写呢？我们公司用的是Abaqus…。
Abaqus通常用SI（米制）书写。同样的1米跌落设置：
** --- 初速度 ---
*INITIAL CONDITIONS, TYPE=VELOCITY
product_set, 3, -4.43
** --- 重力 ---
*DLOAD
product_set, GRAV, 9.81, 0., 0., -1.
** --- 刚体地面（解析刚体面）---
*RIGID BODY, REF NODE=floor_rp,
ANALYTICAL SURFACE=floor_surf
*SURFACE, TYPE=PLANAR, NAME=floor_surf
DATA LINE: 0., 0., 0., 0., 0., 1.
** --- 通用接触（推荐）---
*CONTACT, OP=NEW
*CONTACT INCLUSIONS, ALL EXTERIOR
** --- 分析步（5 ms分析）---
*DYNAMIC, EXPLICIT
, 0.005Abaqus的 General Contact 非常方便，会自动将所有外表面注册为接触候选。这与LS-DYNA的 AUTOMATIC_SINGLE_SURFACE 概念类似。用解析刚体面定义地面，则无需地面的网格，可以提高计算效率。
单元选择与网格策略
应该选择哪种单元类型呢？我在实体单元和壳单元之间犹豫。
跌落分析中根据零件的厚度和变形性质来区分使用：
| 零件 | 推荐单元 | LS-DYNA | Abaqus | 理由 |
|---|---|---|---|---|
| 外壳（1〜3 mm厚） | 壳单元 | ELFORM=2（BT） | S4R | 板厚方向5个积分点确保弯曲精度 |
| 基板（0.8〜1.6 mm） | 壳单元 | ELFORM=16（完全积分） | S4 | 准确捕捉翘曲变形 |
| IC/BGA（焊点连接） | 实体单元 | ELFORM=1（减缩积分） | C3D8R | 需要评估三维应力状态 |
| 橡胶垫圈 | 实体单元 | ELFORM=-1（完全积分） | C3D8H（混合法） | 避免不可压缩橡胶的体积自锁 |
| 发泡缓冲垫 | 实体单元 | ELFORM=1 | C3D8R | 应对大压缩变形 |
网格尺寸需要相对于冲击波的波长足够精细。经验法则是接触面附近单元尺寸为 0.5〜1.0 mm，其他部位为2〜5 mm左右，厚度方向壳单元建议5个积分点，实体单元至少3层。
如果用减缩积分出现沙漏模式怎么办？
なった
詳しく
報告