应力奇异点与网格收敛——有限元法中应力发散原理及实际对策

不是缺陷。弹性力学理论解本身在某些情况下就预测了“应力无限大”。90度的凹角（L形内角270°的角）和裂纹尖端就是典型例子。

对，这正是关键所在。FEM是求弹性理论近似解的，所以网格越细就越接近“理论上的无限大”。也就是说 $\sigma_{\max}$ 在网格尺寸 $h \to 0$ 时会无限增大。通常的网格收敛性检查——用三个级别的网格确认应力是否收敛于一个定值——在此不成立。

正是如此。FEM初学者用最大von Mises应力与许用应力比较——这就是包含奇点的模型中典型的过度设计模式。网格粗的人得到“OK”，网格细的人得到“NG”，结果会因网格负责人的主观判断而改变。这在V&V（验证与确认）上完全是不合格的。

有Williams（1952）推导出的经典特征值分析。设楔形区域顶点距离为 $r$，角度为 $\theta$，则应力场可表示为：

这里 $K$ 是与应力强度相关的系数，$f_{ij}(\theta)$ 是角度相关函数，而 $\lambda$ 是Williams特征值。这个 $\lambda$ 取决于楔形张开角 $2\alpha$。关键在于，当 $\lambda < 1$ 时，$r \to 0$ 则 $\sigma \to \infty$。

求解Williams特征值方程可知，内角270°（凹角）时 $\lambda \approx 0.545$。也就是说应力以 $\sigma \propto r^{-0.455}$ 的量级发散。裂纹尖端（内角360°）时则是众所周知的 $\lambda = 0.5$，$\sigma \propto r^{-1/2}$。

没错。所以在实际工作中，通过添加圆角R（倒圆）来消除角部，从而避免奇异性。但是，即使CAD模型有圆角，如果网格太粗以至于无法充分表现角部形状，数值上仍可能残留奇异性，需要注意。

准确地说，远离奇点位置的应力确实会收敛。问题在于奇点正附近的最大应力。我们来看标准FEM的误差估计式：

这里 $k$ 是单元多项式次数，$\lambda$ 是Williams特征值。对于光滑解，$\lambda \geq k$ 时可获得 $h^k$ 的收敛率。但存在奇点时 $\lambda < 1 < k$，所以收敛率被限制在 $h^\lambda$。即使使用二次单元，也只能达到 $h^{0.545}$ 左右，而不是 $h^2$。

对于全局能量范数误差来说，并非完全没用，只是改善缓慢。但是，奇点处的应力值本身，即使提高单元次数或细化网格也不会收敛——这是本质上的局限。这不是FEM的缺陷，而是物理模型（线性弹性体的锐角角部）所具有的本质特性。

正是如此。所以“用哪里的应力来评价”是奇点问题的本质。看到最大应力的云图就反射性地判断“红色区域危险！”，这在有奇点的模型中是完全错误的。

列举一些代表性的：

• 凹角 — L形支架、T形连接处无圆角的角部。实际工作中遭遇频率最高
• 裂纹尖端 — 断裂力学模型中故意引入的。$\lambda = 0.5$，奇异性最强
• 异种材料界面的端部 — 杨氏模量不同的材料接合端。称为双材料奇异性
• 集中载荷作用点 — 点载荷、线载荷理论上接触面积为零 → 应力无限大
• 约束条件的角部 — 固定端与自由端边界处的角。会产生数值奇异性
• V型缺口 — 焊接焊趾端部、键槽底部等

是的。如果用“单节点集中载荷”来替代螺栓紧固的座面或轴承的接触面，则该节点附近的应力在物理上没有意义。至少应将载荷分布到面上，或将载荷作用点附近的应力排除在评价对象之外，这是铁则。

问得好。断裂力学中评价的不是“应力值”而是“应力强度因子 $K$”或“$J$ 积分”，所以不使用奇点处的应力值本身。但是，为了用FEM高精度求得 $K$，有一种技巧叫做四分之一节点单元。

这是Barsoum（1976）提出的方法。在8节点四边形二次单元中，将裂纹尖端侧的中节点偏移到边的四分之一位置。这样，单元内的位移场会自然地出现 $\sqrt{r}$ 项，从而能在单元层面准确捕捉裂纹尖端的 $r^{-1/2}$ 奇异性。

Abaqus或Ansys中内置了在裂纹尖端周围配置这种四分之一节点单元的功能。结合后面要讲的J积分或围道积分法，就能高精度求得 $K_I, K_{II}, K_{III}$。

hp细化是有效的。向奇点方向按几何级数（geometric grading）细化网格，同时提高单元的多项式次数。根据Babuska等人的研究，通过这种组合可以实现指数收敛：

这里 $N$ 是自由度数。与通常的h细化只能代数收敛（$N^{-\lambda/d}$）相比，hp细化是指数收敛，效率极高。不过实现复杂，支持的求解器有限。StressCheck或基于p-FEM的求解器是代表例子。

主要有三种方法：

1. 位移外推法 — 利用裂纹尖端附近节点位移，根据 $K = \lim_{r \to 0} \sigma\sqrt{2\pi r}$ 的关系进行外推。简单但网格依赖性高
2. J积分法（等效区域积分） — 沿围绕裂纹尖端的路径（围道）计算能量释放率 $G$。$K_I = \sqrt{E' G}$（$E'$ 因平面应力/应变而异）。具有路径无关性，精度高
3. 相互作用积分法（Interaction Integral） — 利用辅助场与实际场的J积分相互作用，单独分离提取 $K_I, K_{II}, K_{III}$。混合模式问题必备

正是如此。所以不使用奇点正上方的应力值，而是用稍远的多个围道的平均值来求 $J$。Abaqus中通过 *CONTOUR INTEGRAL 可以自动计算5~10个围道。越外侧的围道受网格影响越小，越稳定。

ASME锅炉及压力容器规范第VIII卷第2册规定的应力线性化（stress linearization）是行业标准。这是一种将板厚方向的应力分布分解为“膜应力 $\sigma_m$”和“弯曲应力 $\sigma_b$”的方法。

这里 $t$ 是板厚，$x$ 是板厚方向坐标。膜应力是均匀拉伸截面的分量，弯曲应力是正反面符号反转的分量。奇点引起的峰值应力 $\sigma_p$（线性化后残留的非线性分布残差）被排除在一次评价对象之外。

在“一次应力评价”中可以忽略。但疲劳评价中峰值应力很重要。ASME Div.2 Part 5.5中，疲劳分析使用全范围应力（膜+弯曲+峰值）。此时奇点的峰值应力在物理上没有意义，因此切换到带圆角R的子模型来评价实际的应力集中才是正确做法。