应力奇异点与网格收敛 — FEM应力发散的原理与实务对策

不是BUG。弹性力学的理论解本身就预测"应力无限大"的情况。90度的重入角（L字内角270°）或裂纹尖端就是典型。

这正是关键。FEM求的是弹性理论的近似解，网格越细越接近"理论无穷大"。所以 $\sigma_{\max}$ 随着 $h \to 0$ 而无限增长。通常的网格收敛性检查——用3个网格密度验证应力收敛——根本不成立。

完全正确。用最大von Mises应力与允许应力比较——这是包含奇异点模型的典型过度设计陷阱。粗网格通过，细网格不通过，结果完全取决于网格控制者的做法。从V&V角度这是完全错误的。

Williams（1952）推导了经典的特征值分析。在楔形区域，以顶点为原点，距离为 $r$、角度为 $\theta$，应力场表示为：

这里 $K$ 是应力强度系数，$f_{ij}(\theta)$ 是角度函数，$\lambda$ 是Williams特征值，它取决于楔形开口角。关键是当 $\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模型有R，如果网格太粗不足以表现R形状，数值上特异性仍然存在。

准确说，奇异点远处的应力是收敛的。问题在奇异点附近的最大应力。标准FEM误差估计式是：

$k$ 是单元多项式次数，$\lambda$ 是Williams特征值。光滑解时 $\lambda \geq k$ 得到 $h^k$ 收敛率。但有奇异点时 $\lambda < 1 < k$，收敛率被限制到 $h^\lambda$。2次单元也只能得到 $h^{0.545}$ 的收敛，不是 $h^2$。

对全局能量误差有改进但很慢。关键是奇异点处应力值本身不收敛——这是固有限制。不是FEM的问题，而是线性弹性模型在尖锐角处的本质性质。

完全正确。所以评估应力的位置至关重要。看应力云图的红区就是最大应力，这在含奇异点的模型上是陷阱。

常见的有：

是的。螺栓座面或轴承接触面用"1个节点的集中荷载"代替，那个节点附近的应力是物理无意义的。至少要把荷载分散在一个面，或者不评估荷载加载点附近的应力。

破坏力学不评估应力值，改用应力强度因子 $K$ 或 $J$ 积分。为了精确求 $K$，Barsoum（1976）发明了四分点单元。

8节点二次四边形单元，把裂纹尖端附近的中间节点从边中点移到 $1/4$ 处。这样单元内的位移自动包含 $\sqrt{r}$ 项，能精确捕捉 $r^{-1/2}$ 奇异性。

Abaqus、Ansys都有裂纹尖端附近自动配置四分点单元的功能。配合J积分或轮廓积分法就能高精度求 $K_I, K_{II}, K_{III}$。

有，**hp-细分化**很有效。朝奇异点以几何级数细分网格，同时增加多项式次数。Babuska等的理论证明能达到指数收敛：

普通h细分只有代数收敛 $N^{-\lambda/d}$，hp-细分是指数级。但实现复杂，支持hp的商业求解器很少。StressCheck和p-FEM型求解器能做。

主要3种方法：

对。所以不用在尖端直接算，而是用离尖端稍远的多条路径平均。Abaqus的 *CONTOUR INTEGRAL 自动计算5-10条路径。外层路径受网格影响小，结果更稳定。

ASME Boiler & Pressure Vessel Code Section VIII Division 2 规定的应力线性化是业界标准。把板厚方向的应力分解成膜应力 $\sigma_m$ 和弯曲应力 $\sigma_b$。

$t$ 是板厚，$x$ 是厚度坐标。膜应力是均匀拉伸，弯曲应力上下符号相反。奇异点的峰值应力 $\sigma_p$ 在一次评估时被排除。

一次评估不考虑。但疲劳评估时峰值很重要。ASME Div.2 Part 5.5疲劳分析用 full-range stress（膜+弯曲+峰值）。此时奇异点的峰值物理无意义，要用含圆角R的子模型重算实际的应力集中。

完美理解！加上第0条就是设计阶段就加圆角，根本上避免奇异点。

切换到破坏力学参数评估。主要3种方法：

完全正确。应力发散成 $r^{-1/2}$，但J积分路径无关且有限，$K$ 也唯一确定。从点值转向整体参数是奇异点问题的核心解决方案。

子模型有2个用途：

注意！子模型的截断边界必须远离奇异点。全局模型在那个位置的应力必须已经收敛，否则边界条件不准。

加上第0条：设计时加圆角，根除奇异点！

比较一下主要的：

是。WorkBench指定板厚方向路径，自动算 $P_m, P_b, P_m+P_b, P_m+P_b+Q$，压力容器设计审查时很省事。Abaqus要用 *PATH 定义再Python脚本计算，不如一键便利。

开源方面，CalculiX支持J积分。Code_Aster（法国EDF开发）功能很强，POST_K1_K2_K3命令能抽应力强度因子，还支持XFEM，破坏力学方面不逊于商用。

传统FEM裂纹必须沿单元边界。XFEM在现有网格基础上加入富集函数，不需单元边界对齐裂纹，甚至网格可以不动。

$H(\mathbf{x})$ 表示裂纹面的不连续，$F_\ell$ 是 $\sqrt{r}$ 特异性函数族。形状函数层面就把尖端特异性吸收了，粗网格也能求出高精度的 $K$ 值。

对。这是XFEM的革命性优势。缺点是3D复杂形状、裂纹分叉合并等还有工程化的困难。Abaqus的3D XFEM基本能用，但裂纹扩展限于疲劳Paris律，脆性急速扩展在3D的稳定性仍有问题。

相场法完全不同。不用鋭利裂纹，而用连续的损伤变量 $d \in [0, 1]$ 表示。$d=0$ 健全，$d=1$ 完全破坏。Bourdin等（2000）基于Griffith能量平衡推导。

第一项是劣化的弹性能，第二项是裂纹表面能。$G_c$ 是临界能量释放率，$\ell$ 正则化长度。裂纹的发生、扩展、分叉自动出现，不用预设裂纹路径。

目前（2026年）商用求解器XFEM先行。相场法还主要在研究阶段，虽然Abaqus 2024开始加强了UEL/UMAT的相场支持，但远未成为标配。不过学术论文近5年爆炸增长，5-10年内应该标配化。

见过的失败案例5个：