惩罚法接触定式化
惩罚法接触定式化的理论基础
接触问题的基础
老师，有限元的「接触问题」为什么这么难？
接触问题是本质上非线性的。有三个理由：
1. 状态变化 — 接触面是「接触/非接触」的二值。荷重改变状态
2. 不等式约束 — 不能发生贯通（penetration）：$g \geq 0$
3. 摩擦 — 库仑摩擦包含「粘着/滑移」的状态变化
不等式约束…通常有限元的等式约束（SPC等）在本质上不同，是吧。
通常有限元求解 $[K]\{u\} = \{F\}$ 的等式，但接触需要满足「$g \geq 0$ 且 $g \cdot p = 0$」的KKT条件（Karush-Kuhn-Tucker条件）。$g$ 是间隙，$p$ 是接触压。
惩罚法的原理
惩罚法是最简单的接触处理方法。「贯通的反力与贯通量成正比」：
$k_p$ 是惩罚刚度（接触刚度）。$g_n$ 是法向贯通量。
接触面上放个弹簧的感觉吗？
完全正确。惩罚法是在接触面上配置一个非常硬的弹簧。发生贯通时，反力产生并推回贯通。
惩罚刚度的设置
惩罚刚度 $k_p$ 的设置影响结果：
- $k_p$ 太大 → 条件数恶化。收敛困难
- $k_p$ 太小 → 贯通量过大。精度不足
经验值：$k_p \approx 10 \sim 100 \times E \cdot A / L$（接触面刚度的量级）。大多数求解器自动计算。
交给自动计算安全吗？
大多数情况下，是的。Abaqus和Ansys的自动惩罚刚度是个好起点。只有出现问题时才手动调整。
惩罚法的优缺点
| 优点 | 缺点 |
|---|---|
| 实现简单 | 贯通无法完全消除 |
| 不需要增加自由度 | 惩罚刚度设置影响结果 |
| 与显式求解相容性好 | 大的$k_p$会恶化条件数 |
| 大多数求解器默认方法 | — |
「贯通无法完全消除」是最大的弱点吗？
由于 $k_p$ 有限，微小贯通无法避免。只要贯通量低于板厚的1%，工程上就可接受。超过这个值就考虑用拉格朗日乘数法。
小结
要点：
- 接触是不等式约束的非线性问题 — 状态变化＋摩擦
- 惩罚法 = 接触面上的弹簧 — $F = k_p \cdot g_n$
- $k_p$ 设置是关键 — 太大收敛困难，太小贯通过大
- 大多数求解器自动计算 — 只在有问题时手动调整
- 贯通无法完全消除 — 目标是板厚的1%以下
Courant 1943年的罚则法
惩罚法的数学起源来自Richard Courant在1943年发表的论文「Variational methods for the solution of problems of equilibrium and vibrations」。在约束变分问题上乘以大系数（惩罚），将其转化为无约束问题的想法，最初用于椭圆型偏微分方程的边界条件处理。1970年代Oden & Kim(1977)等人系统地应用于接触有限元。
惩罚法接触定式化的数值计算方法
惩罚法的实现
在各求解器中如何设置惩罚法？
Abaqus
```
*CONTACT PAIR, INTERACTION=contact_prop
slave_surface, master_surface
*SURFACE INTERACTION, NAME=contact_prop
*SURFACE BEHAVIOR, PENALTY
```
Abaqus默认是惩罚法。用 *SURFACE BEHAVIOR, PENALTY 明确指定。
Nastran
```
BCTABLE, ...
BCPARA, CTYPE, UGLYPEN $ 惩罚法
```
Nastran的SOL 101/106/400中定义接触。
Ansys
```
MP, MU, cid, 0.3 ! 摩擦系数
KEYOPT, cid, 2, 1 ! 惩罚法
```
或在Workbench接触设置中设置Formulation=Penalty。
LS-DYNA
```
*CONTACT_AUTOMATIC_SURFACE_TO_SURFACE
$ 惩罚法为默认
```
LS-DYNA的所有接触都基于惩罚法。显式方法最自然。
接触检测
惩罚法的步骤（每次迭代）：
1. 接触检测 — 在主面上找从节点的最近点
2. 间隙计算 — $g_n$ = 从节点与主面的距离
3. 贯通判定 — $g_n < 0$ 则接触。反力 $F_n = k_p |g_n|$
4. 摩擦力计算 — $F_t = \mu F_n$（滑移）或粘着
5. 反映到全局方程 — 惩罚力加到右边
接触检测占了大部分计算时间吗？
在大规模模型（数十万接触单元对）中，接触检测占计算时间的30～50%。使用Bucket sort或KD树等高速搜索算法。
小结
惩罚系数的选择准则
接触惩罚系数ε的选择是惩罚法实用上的最重要课题。ε太小浸透量超过容许值，太大刚性矩阵条件数恶化收敛变慢。LS-DYNA默认设置采用「SOFT=0」算法，自动将接触面最小单元刚度k_e（=E×t×A）的0.1倍设为ε，异种材料接触（如钢-树脂）时推荐手动调整。
惩罚法接触定式化的实务应用
接触分析的实务
请教接触分析实务的要点。
主面和从面的选择
为什么软的一侧要设为从面？
惩罚法中从节点不能穿过主面。反过来就是主面会「扎进」软的从面。将硬面作为主面，软面作为从面，可以避免「扎穿」。
接触面的网格
收敛的技巧
改善接触分析收敛的方法：
1. 初始间隙设为零 — 从一开始就接触的状态开始
2. 分阶段施加荷重 — 初始增量要小
3. 接触稳定化 — Abaqus的*CONTACT STABILIZATION
4. 调整惩罚刚度 — 从自动改为手动
5. 分步引入摩擦 — 先$\mu=0$建立接触，再加摩擦
实务检查清单
「验证贯通量」是惩罚法最重要的检查，对吧。
贯通过大则接触压精度降低。整个接触面的贯通量用后处理器检查。
摩托车头盔冲击分析
摩托车安全标准ECE 22.06（2020年版）规定，可用CAE评估斜面冲击时的脑加速度（HIC值）。EPS（泡沫塑料）衬垫与纤维增强外壳的接触广泛使用惩罚法。Shoei公司的公开数据表明，使用ABAQUS Explicit的惩罚接触进行HIC分析，精度达到实测值±12%以内，已被2022年后的新品认证采用。
惩罚法接触定式化的软件比较
接触求解器比较
Ansys推荐增强拉格朗日吗？
Ansys默认是增强拉格朗日。惩罚法+迭代修正贯通，纯惩罚法的贯通更小。Abaqus默认也是惩罚法，但KINEMATIC（拉格朗日乘数）也广泛使用。
选择指南
LS-DYNA惩罚法的历史
John Hallquist在1976年在LLNL开发的DYNA2D中的接触实现最初只是简单的节点-表面型惩罚法。1989年的LS-DYNA3D V900版本增加了segment-to-segment接触，能更准确地评估面压分布。最新的R15版本允许单个分析模型定义最多10,000个接触，远超汽车全车碰撞分析的需求。
惩罚法接触定式化的先端研究
Mortar法
Mortar法是克服惩罚法弱点（网格依赖性、贯通）的最新接触方法。用弱形式（积分形式）施加接触条件，对网格不匹配有强大容错能力。Abaqus的SURFACE TO SURFACE接触就是基于Mortar。
IGA接触
等几何分析（IGA）用NURBS基函数描述接触面。曲面接触光滑，不会出现普通有限元网格依赖的接触压振荡。
摩擦的微观模型
对表面微观粗糙度建模来「预测」摩擦系数的研究。多尺度接触力学。
小结
GPU惩罚接触的高速化
2020年以后，将惩罚接触的浸透检测和力计算用GPU并行执行的架构进入实用阶段。Ansys LS-DYNA GPU版（2022年发布）在CUDA内核中实现了bucket sort近邻搜索，在100万单元规模的汽车碰撞分析中达到CPU版最大8倍的加速。GPU接触与CPU版精度差异小于0.1%，官方保证。
惩罚法接触定式化的故障排除
接触分析的故障
接触分析常见的故障有哪些？
不收敛
接触分析不收敛是最常见的故障。对策：
1. 初始增量减小（0.1→0.01→0.001）
2. 使用接触稳定化（*CONTACT STABILIZATION）
3. 清除初始间隙（从已接触状态开始）
4. 分步引入摩擦（$\mu=0$→$\mu=0.1$→$\mu=0.3$）
5. 降低惩罚刚度（贯通增加但更容易收敛）
贯通量过大
惩罚刚度太小。提高 $k_p$，或切换到增强拉格朗日法。
接触压呈棋盘图案
一阶单元（TET4、HEX8完全积分）的接触面压力振荡。改用C3D10M（Abaqus），或二阶单元降积分。
从节点穿入主面
主/从面选择反了。检查「硬面=主面」的原则。
小结
能量误差和沙漏
惩罚接触分析的典型失败是「hourglass（沙漏）能量过大」。低阶降积分单元与惩罚接触结合时，接触力会激发沙漏变形，导致数值发散。2008年Livermore Software的基准报告表明，当ERODED_SINGLE_SURFACE接触的沙漏能量超过总变形能的15%时，结果信度严重下降，需要加强沙漏控制或改变接触定式。
相关主题
价值
更详细
错误