中央差分法（陽解法）
理论与物理
显式解法是什么
不需要解联立方程？那看起来很快。
单步计算速度快几个数量级。但有稳定条件限制，时间步长必须非常小（$\Delta t \propto L_{min} / c$，$c$ 是波速）。最适合冲击这类短时间现象。
中心差分法的算法
对于运动方程 $[M]\{\ddot{u}\} = \{F\} - [K]\{u\} - [C]\{\dot{u}\}$：
1. 加速度计算: $\{\ddot{u}_n\} = [M]^{-1}(\{F_n\} - \{F_{int,n}\})$
2. 速度更新: $\{\dot{u}_{n+1/2}\} = \{\dot{u}_{n-1/2}\} + \Delta t_n \{\ddot{u}_n\}$
3. 位移更新: $\{u_{n+1}\} = \{u_n\} + \Delta t_{n+1/2} \{\dot{u}_{n+1/2}\}$
如果 $[M]^{-1}$ 是对角矩阵，那么求逆就只是取每个对角元素的倒数…计算超快！
所以显式解法必须使用集中质量矩阵（对角）。如果使用一致质量矩阵（非对角），求逆计算成本高，会失去显式解法的优势。
稳定条件（CFL条件）
显式解法的稳定时间增量：
$$ \Delta t \leq \frac{L_{min}}{c} = \frac{L_{min}}{\sqrt{E/\rho}} $$
显式解法的稳定时间增量：
$L_{min}$ 是最小单元尺寸，$c$ 是材料的波速（弹性波速度）。
钢材（$c \approx 5000$ m/s）且 $L_{min} = 5$ mm 时，$\Delta t \leq 1 \mu s$…非常小！
所以分析1秒需要超过100万步。冲击（几毫秒）只需几千步，但准静态（几秒）则步数巨大。显式解法适用于短时间的动态现象。
显式解法 vs. 隐式解法
| 特性 | 显式解法 | 隐式解法 |
|---|---|---|
| 单步计算 | 非常轻量 | 重（联立方程） |
| 时间步长 | 非常小（$\mu s$量级） | 大（$ms \sim s$） |
| 稳定条件 | 有（CFL条件） | 无（无条件稳定） |
| 适用场景 | 冲击、爆炸、碰撞 | 准静态、振动 |
| 接触 | 容易 | 困难 |
| 非线性 | 自然处理 | 需要迭代 |
接触容易处理是个很大的优点呢。
显式解法中，接触检测和力的计算在每个时间步独立进行，没有收敛问题。对于像汽车碰撞这样同时发生大量接触的问题，显式解法具有压倒性优势。
总结
我来整理一下中心差分法（显式解法）。
要点：
- 不解联立方程 — 若 $[M]^{-1}$ 为对角阵则只需取倒数
- $\Delta t \leq L_{min} / c$ — CFL条件。由最小单元决定
- 最适合短时间动态现象 — 冲击、爆炸、碰撞
- 接触容易处理 — 无收敛问题
- 必须使用集中质量矩阵 — 非对角矩阵则无法提速
显式解法是“针对短时间剧烈现象”的算法呢。
LS-DYNA的全部分析都是显式解法。Abaqus/Explicit也是显式解法。汽车的碰撞安全没有这个方法就无法实现。
中心差分法诞生于1695年的微积分
莱布尼茨在1695年发表的有限差分概念被应用于动力学，产生了中心差分法，它用 a(t)＝[u(t+Δt)−2u(t)＋u(t−Δt)]/Δt² 来近似加速度。这个二阶精度格式从冯·诺依曼稳定性分析导出了临界时间步长 Δtcr＝2/ωmax，约束其不超过最小单元的固有频率。该算法于1970年代在LS-DYNA的前身代码DYNA3D中首次大规模实现于FEM。
各项的物理意义
- 惯性项（质量项）：$\rho \ddot{u}$，即“质量×加速度”。您有过急刹车时身体被向前甩出去的经历吗？那种“被带走的感觉”正是惯性力。物体越重越难启动，一旦动起来也越难停止。地震时建筑物摇晃，也是因为地面突然移动，而建筑物的质量“被落下”。静力分析中此项设为零，那是假设“因为缓慢施力所以加速度可以忽略”。对于冲击载荷或振动问题，此项绝对不能省略。
- 刚度项（弹性恢复力）：$Ku$ 或 $\nabla \cdot \sigma$。拉弹簧时会感觉到“想恢复原状的力”吧？那就是胡克定律 $F=kx$，也是刚度项的本质。那么提问——铁棒和橡皮筋，用相同的力拉，哪个伸得更长？当然是橡皮筋。这种“难拉伸程度”就是杨氏模量 $E$，它决定了刚度。常见的误解：“刚度高＝强度高”是不对的。刚度是“难变形程度”，强度是“难破坏程度”，是不同的概念。
- 外力项（载荷项）：体积力 $f_b$（重力等）和表面力 $f_s$（压力、接触力等）。可以这样想——桥上卡车的重量是“作用在整个内部上的力”（体积力），轮胎压路面的力是“只作用在表面上的力”（表面力）。风压、水压、螺栓紧固力…全都是外力。这里容易犯的错误：弄错载荷方向。本想“拉伸”却成了“压缩”——听起来像笑话，但在3D空间中坐标系发生旋转时确实会发生。
- 阻尼项：瑞利阻尼 $C\dot{u} = (\alpha M + \beta K)\dot{u}$。弹一下吉他弦试试。声音会一直响吗？不，会逐渐变小对吧。因为振动能量通过空气阻力和弦的内部摩擦变成了热。汽车的减震器也是同样原理——故意吸收振动能量来改善乘坐舒适性。如果阻尼为零会怎样？建筑物在地震后会一直摇晃不停。实际上不会这样，所以设置适当的阻尼很重要。
假设条件与适用范围
量纲分析与单位制
| 变量 | SI单位 | 注意事项·换算备忘 |
|---|---|---|
| 位移 $u$ | m（米） | 以mm输入时，载荷·弹性模量也需统一为MPa/N系 |
| 应力 $\sigma$ | Pa（帕斯卡）= N/m² | MPa = 10⁶ Pa。与屈服应力比较时注意单位制不一致 |
| 应变 $\varepsilon$ | 无量纲（m/m） | 注意工程应变与对数应变的区别（大变形时） |
| 弹性模量 $E$ | Pa | 钢：约210 GPa，铝：约70 GPa。注意温度依赖性 |
| 密度 $\rho$ | kg/m³ | mm制时为tonne/mm³（钢为 = 10⁻⁹ tonne/mm³） |
| 力 $F$ | N（牛顿） | mm制用N，m制也用N统一 |
数值解法与实现
显式解法的实现
显式解法的FEM代码是如何实现的？
主循环极其简单：
```
for each time step:
1. 内力计算: F_int = assemble(sigma B V) $ 单元循环
2. 外力计算: F_ext = boundary_conditions()
3. 加速度: a = (F_ext - F_int) / M_diag $ 对角除法
4. 速度更新: v = v + dt * a
5. 位移更新: u = u + dt * v
6. 接触检查: contact_detection_and_force()
7. 时间增量更新: dt = CFL * L_min / c
```
第3步的“对角除法”是关键呢。
是的。这里是显式解法速度的源泉。隐式解法中第3步是“求解 $n \times n$ 的联立方程”，计算量差几个数量级。
稳定时间增量的管理
稳定时间增量由所有单元中最小的值决定：
只要有一个单元很小，整体的 $\Delta t$ 就会变小。
因为一个坏单元导致整体变慢？
网格质量直接关系到计算速度。长宽比差的单元或压扁的单元 $L_e$ 非常小，会减小 $\Delta t$。显式解法中网格质量是计算效率的关键。
质量缩放
人为增加小单元的质量以增大 $\Delta t$ 的质量缩放：
增加 $\rho$ 会降低 $c$，从而增大 $\Delta t$。是准静态问题显式解法中必备的技术。
增加质量不会改变惯性效应吗？
会改变。所以需要确认动能/内能 < 5%。满足这个条件就可以视为“准静态”。
求解器
相关主题この記事の評価 ご回答ありがとうございます！ 参考に なった もっと 詳しく 誤りを 報告 |
|---|