运动学基础 — 位移、速度、加速度与有限元时间积分
目录
1. 运动学在CAE中的角色
跌落试验仿真的时候，软件需要设置初始速度，这个怎么算？产品从1米高度自由落下……
直接用自由落体公式：$v = \sqrt{2gh}$。1米高度的话，$v = \sqrt{2 \times 9.81 \times 1} \approx 4.43$ m/s，方向向下。在LS-DYNA或Abaqus里设初始条件，给整个模型施加这个速度就行。
运动学（Kinematics）研究物体运动的几何描述，不涉及力的原因。在CAE中，运动学关系贯穿于：
- 初始条件设置：跌落试验、碰撞分析的初速度
- 时间积分：从加速度推算速度和位移的数值方法
- 应变率计算：$\dot{\varepsilon} = \frac{d\varepsilon}{dt}$，影响率相关材料（金属高速变形）
- 运动边界条件：强制位移输入（地震激励）
2. 位移：运动的基本描述
位移 $\vec{u}$ 是物体从初始位置到当前位置的向量：
在FEM中，位移是最基本的未知量。对于三维问题，每个节点有3个位移自由度：$u_x, u_y, u_z$（对于固体单元）；梁壳单元还有3个转角自由度：$\theta_x, \theta_y, \theta_z$。
位移场在单元内通过形函数插值：
这就是FEM"有限元"名称的由来：用有限个节点的位移，通过形函数近似描述连续的位移场。
3. 速度：位移对时间的导数
速度是位移对时间的一阶导数：
速度在CAE中的主要作用：
- 动能计算：$T = \frac{1}{2}\int_V \rho |\dot{u}|^2 \, dV$
- 阻尼力：$F_{\text{damp}} = [C]\{\dot{u}\}$（与速度成正比）
- 初始速度条件：显式动力学仿真的起始状态
- 应变率：$\dot{\varepsilon}_{ij} = \frac{1}{2}(v_{i,j} + v_{j,i})$
应变率这个东西在哪种工程问题里特别重要？平时做静力分析用不到吧？
静力分析确实可以忽略。但高速碰撞就必须考虑了。金属在高应变率下的屈服强度会显著升高——这叫"应变率强化效应"。比如低碳钢在应变率1000/s时的屈服强度可能是准静态的1.5到2倍。汽车碰撞吸能盒的材料模型（如Johnson-Cook本构）里就有应变率项，不考虑的话吸能预测会偏差很大。
4. 加速度：速度对时间的导数
加速度是速度对时间的一阶导数（位移的二阶导数）：
加速度在工程分析中的意义：
| 应用场景 | 典型加速度量级 | CAE分析类型 |
|---|---|---|
| 汽车碰撞乘员舱 | 30-50 g | 显式动力学（LS-DYNA） |
| 电子产品跌落 | 500-2000 g | 显式动力学（Abaqus Explicit） |
| 航空发动机振动 | 5-20 g | 模态/谐响应（Nastran） |
| 地震响应 | 0.1-1 g（PGA） | 瞬态/响应谱（ANSYS） |
| 火箭发射 | 3-6 g（轴向） | 随机振动（Nastran） |
加速度信号（g-t曲线，即加速度-时间历程）是碰撞/振动试验的核心测量量，CAE仿真结果必须与之对标。
5. 匀变速运动方程组
当加速度为常数时（$a = \text{const}$），运动学方程组为：
这三个方程组成了分析匀变速运动的完整工具集。在实际CAE仿真中，这组方程用于：
- 估算跌落试验的冲击速度（$v = \sqrt{2gh}$）
- 估算制动距离（$x = v_0^2 / (2a)$）
- 验证仿真结果的量级合理性（"Back-of-envelope" 检验）
实例计算：一款智能手机从办公桌（高度 $h = 0.75$ m）跌落到地板。冲击瞬间的速度：
这个速度就是Abaqus Explicit跌落分析的初始速度边界条件。
6. 工程实例：跌落试验初速度计算
笔记本电脑的跌落标准是什么？仿真时怎么设置工况？
MIL-STD-810和IEC 60068-2-31是常用标准。消费电子通常测试0.9m到1.2m高度，跌落方向包括6面、12棱、8角——共26种工况。仿真时对应地要建26个分析任务。不过一般先做最严苛的角跌（Corner drop），因为那时候冲击能量集中在最小接触面积上，应力最大。
| 跌落高度 | 冲击速度 | 等效减速度（接触时间10ms） | 典型适用标准 |
|---|---|---|---|
| 0.5 m | 3.13 m/s | ~32 g | 消费电子（低端） |
| 0.9 m | 4.20 m/s | ~43 g | 笔记本电脑（IEC 60068） |
| 1.2 m | 4.85 m/s | ~49 g | 工业手持设备（MIL-STD-810） |
| 1.5 m | 5.42 m/s | ~55 g | 军用设备 |
| 2.0 m | 6.26 m/s | ~64 g | 特种设备 |
Abaqus跌落分析典型设置：
地板建模为刚性面（Rigid surface），接触属性设置摩擦系数 $\mu = 0.3$（典型橡胶-钢）。分析时间通常取 1~5 ms（覆盖整个冲击脉冲）。
7. FEM时间积分方法
动力学问题的核心挑战是：如何从已知的 $\{u\}_n, \{\dot{u}\}_n, \{\ddot{u}\}_n$ 推进到下一时刻 $\{u\}_{n+1}$？这是时间积分问题。
Newmark-β 法（隐式）：
常用参数：$\beta = 0.25, \gamma = 0.5$（平均加速度法，无条件稳定）。
中心差分法（显式）：
显式法不需要求解线性方程组，但时间步长受稳定条件限制（条件稳定）。
8. 显式与隐式积分的选择
碰撞用显式，振动用隐式——这个规律是怎么来的？有没有反过来用的情况？
规律来自效率考量。显式每步计算很便宜，但步长极小（微秒级），适合碰撞这种分析时间短（100ms量级）的问题。振动分析时间可能长达几秒甚至几十秒，用显式要算几千万步，代价太高——这时隐式法步长可以取得很大（毫秒级），总步数少得多。反过来用的例子：焊接分析用显式（局部高速热源移动）；航空结构鸟撞有时也从显式切换到隐式来分析后期低速大变形阶段。
| 特性 | 显式积分（中心差分） | 隐式积分（Newmark-β） |
|---|---|---|
| 稳定性 | 条件稳定（需满足 $\Delta t < \Delta t_{crit}$） | 无条件稳定（$\beta \geq 0.25$） |
| 每步计算量 | 小（无需求解方程组） | 大（需解 $[K_{\text{eff}}]\{u\}=\{F_{\text{eff}}$}） |
| 时间步长 | $10^{-6}$ ~ $10^{-5}$ s | $10^{-4}$ ~ $10^{-2}$ s |
| 适用问题 | 碰撞、爆炸、跌落（短时动态） | 地震、振动、疲劳（长时动态） |
| 代表软件 | LS-DYNA, Abaqus Explicit | Nastran SOL112, Abaqus Standard |
9. 数值稳定性与时间步长
显式积分的临界时间步长（CFL条件）：
其中 $L_e$ 是网格中最小单元的特征长度，$c_s$ 是材料中的纵波速度。
常见材料的纵波速度：
| 材料 | 弹性模量 E (GPa) | 密度 ρ (kg/m³) | 纵波速度 $c_s$ (m/s) |
|---|---|---|---|
| 钢（结构） | 210 | 7850 | 5172 |
| 铝合金 | 70 | 2700 | 5092 |
| 高强度铝 | 73 | 2800 | 5109 |
| CFRP（纵向） | 130 | 1600 | 9014 |
| 橡胶（软） | 0.001 | 1100 | 30 |
实践中，LS-DYNA会自动计算所有单元的临界时间步，取最小值，再乘以安全系数（默认0.9）：
这就解释了为什么网格中存在一个极小单元（如楔形单元、锐角单元）会严重拖慢整个分析——这个"时间步杀手"会把整体步长拉低几十倍。网格质量检查（最小单元长度）是碰撞仿真前处理的重要步骤。
如果确实有很小的单元，但又不想重新划网格，有没有办法提高时间步？
有两种技术：一是"质量缩放（Mass Scaling）"——人为增加小单元的密度，让 $c_s$ 降低，从而提高临界时间步。但这会引入虚假惯性力，只在小单元占总质量比例很小时才可以接受。二是子循环（Sub-cycling）——小单元用小步长，大部分单元用大步长，在特定软件（如Abaqus Explicit）中支持，实现代价较高。
总结
运动学三量（位移、速度、加速度）的微积分关系是FEM时间积分的数学基础。从自由落体公式到中心差分积分格式，本质上都是同一套微分关系的应用。工程师需要根据问题的时间尺度（毫秒vs秒）选择显式或隐式积分，并通过网格质量控制保证数值稳定性。