中央差分法（陽解法）
中央差分法（陽解法）的理论基础
陽解法是什么
不求解联立方程组？这听起来很快。
每步的计算快得多。但有稳定条件，时间步长极小（$\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条件）
陽解法的稳定时间增量：
$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²。这个2阶精度格式通过冯·诺依曼稳定性分析得出临界时间步Δtcr＝2/ωmax，因此受最小单元固有频率的限制。1970年代LS-DYNA的前身代码DYNA3D首次将该算法大规模用于FEM。
中央差分法（陽解法）的数值计算方法
陽解法的实现
陽解法的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%需要被确认。如果满足这个条件，就可以认为是「准静态」。
求解器
| 求解器 | 陽解法 |
|---|---|
| LS-DYNA | 陽解法是主力。碰撞安全的行业标准 |
| Abaqus/Explicit | Abaqus Standard的Explicit版本 |
| Ansys LS-DYNA | Ansys集成的LS-DYNA |
| PAM-CRASH | ESI的碰撞分析求解器 |
| RADIOSS | Altair的陽解法求解器 |
总结
陽解法的实现，总结一下。
要点：
- 主循环极其简单 — 对角矩阵逆数用于加速度计算
- 最小单元支配 $\Delta t$ — 网格质量是计算速度的关键
- 质量缩放可适用于准静态 — 运动能量 < 5%需验证
- LS-DYNA是碰撞安全的行业标准 — 所有汽车制造商使用
- Abaqus/Explicit是Abaqus家族的陽解法 — Standard→Explicit切换
质量缩放使计算快10倍
陽解法的时间步由最小单元尺寸支配，少数微小单元会使整体计算成本提高100倍以上。LS-DYNA的质量缩放（MASS SCALING）通过给微小单元加仓虚质量来提高时间步，当追加质量小于总质量的1-2%时，对解的影响在容许范围内。但在瞬间碰撞等惯性力支配的问题中，3%超会出现明显误差。
中央差分法（陽解法）的实务应用
陽解法的实务应用
陽解法在实务中如何应用？
能量平衡的确认
陽解法结果妥当性确认中最重要的是能量平衡：
检查项目：
- 全能量 — 保存（数%以内变动）
- 沙漏能量 — 全能量的5%以下
- 运动能量（准静态时） — 全能量的5%以下
- 接触能量 — 不为负（穿透迹象）
能量平衡匹配就能信任结果吗？
是必要条件但不充分。能量保存也会出现非物理的变形形式。一定要可视化变形。
实务检查清单
「能量平衡确认」是陽解法最关键的检查。
能量「不会说谎」。有问题一定会在能量平衡中显现。所有陽解法分析都要确认能量输出。
工厂板金压力中陽解法也发挥作用
汽车发动机盖的压力成形分析中，通过将模具速度从实际0.5m/s加速到分析用的5m/s的「工程加速法」与质量缩放结合，用Autoform等增量陽解法代码实现了单工程15〜30分钟计算。但加速比超过10倍时，动态效应会导致回弹量与实测偏离，所以奥迪、宝马等规定内部基准为加速比5倍。
中央差分法（陽解法）的软件比较
陽解法的工具
请比较陽解法求解器。
LS-DYNA压倒性优势吗？
汽车碰撞安全中LS-DYNA是全球所有OEM/Tier1的标准。为EuroNCAP、FMVSS等规格开发的模型（WorldSID、THOR等人体模型）都是LS-DYNA格式提供。
选择指南
「碰撞安全 = LS-DYNA」是无法改变的。
LS-DYNA的30年以上实績和人体模型积累是其他公司无法复制的。碰撞安全新参与者必须与现有LS-DYNA基础兼容。
陽解法代码的起源都在LLNL
现代主要陽解法代码（LS-DYNAPAM-CRASHRadioss）都源自劳伦斯利佛莫尔国家实验室的DYNA3D（1976年）。Radioss是Regie在1980年代针对欧洲汽车开发的独立扩展版本，后来在2006年被Altair收购。PAM-CRASH是ESI Group在1980年代针对汽车优化开发的。因为起源相同，基本算法通用，但在接触处理和材料模型实现上各公司各有差别。
中央差分法（陽解法）的前沿研究
陽解法的前沿研究
请讲陽解法的最前沿。
GPU加速
陽解法各单元计算独立（大规模并行可能），与GPU加速相性极好。LS-DYNA GPU版对特定单元类型（HEX8, SHELL）报告CPU的10〜100倍加速。
100倍！碰撞分析从1小时降到1分钟…。
实际上I/O和CPU-GPU数据传输是瓶颈，实效5〜30倍左右。仍然能大幅缩短设计周期。
陽-陰耦合
结构一部分用陽解法（冲击部），其他部分用陰解法（准静态），同时计算的陽-陰耦合。最大化计算效率。LS-DYNA和Abaqus的一部分支持。
IGA（等几何分析）+ 陽解法
NURBS/T-样条基础IGA单元在陽解法中的应用研究。LS-DYNA已实现IGA壳体。CAD形状精确表示和无网格优点。
总结
子步长循环高效处理异尺度问题
流体-结构耦合或异种材料接触问题中，结构和流体的适当时间步可能相差数十倍。子步长循环（或多时间步）技术在粗网格域用大步长、细网格域用小步长计算，并在边界交换信息。LS-DYNA R12以后的ALE-结构耦合分析中该功能改进，官方文档报告计算时间最多可缩短60%。
中央差分法（陽解法）的故障排除
陽解法的故障
请教陽解法的常见故障。
能量未保存
全能量随时间增加。
接触穿透产生能量。对策：
- 提高接触罚函数刚度
- 细化接触面网格
- 改变接触算法（NTS→MORTAR等）
沙漏模式
低减积分单元（HEX8R, SHELL Q4R）中沙漏模式激励。
对策：
- 监视沙漏能量（< 5%）
- 改变沙漏控制类型（粘性型→刚性型）
- 切换到完全积分单元（计算成本增）
时间增量极小
一个小单元限制整体的 $\Delta t$。
对策：
- 改善网格质量（消除/扩大小单元）
- 质量缩放（针对目标单元）
- 合并或删除问题单元
计算失败（负体积）
单元过度变形，体积变零或负。
对策：
- 单元删除（*MAT_ADD_EROSION等）删除过变形单元
- 细化网格分散变形
- 加入材料模型的破坏准则
总结
陽解法故障处理，总结一下。
持续监视能量。这就是陽解法的铁律。
能量是「宇宙的账本」。账本对不上就有问题。在陽解法中，所有问题都反映在能量平衡中。
负能量是负刚度单元的警告
陽解法分析中内部能量突然变负，通常是破坏/删除单元的残余刚度，或接触罚函数导致的过大应变所致。LS-DYNA的*CONTROL_ENERGY中启用SLIDEOPT＝2（接触能量追踪），在RCFORC输出中找出哪个接触对的能量收支崩溃，这是定式的调试步骤。
更详细
错误