Stress Singularities and Mesh Convergence — Principles of Stress Divergence in FEM and Practical Countermeasures

It's not a bug. There are cases where the theoretical solution of elasticity itself predicts "infinite stress". A 90-degree reentrant corner (a 270° internal angle L-shaped corner) or a crack tip are typical examples.

Yes, that's exactly the point. FEM is finding an approximate solution to the elasticity theory, so as the mesh is refined, it approaches the "theoretical infinity". In other words, $\sigma_{\max}$ increases without bound as the mesh size $h \to 0$. The usual mesh convergence check—verifying if the stress converges to a constant value with three mesh levels—does not hold.

Exactly right. A FEM beginner comparing the maximum von Mises stress with the allowable stress—this is the typical over-design pattern for models containing singularities. The person with a coarse mesh gets "OK", the person with a fine mesh gets "NG", so the result changes depending on the mesh operator's discretion. This is completely unacceptable from a V&V standpoint.

There is a classical eigenvalue analysis derived by Williams (1952). Let the distance from the vertex of a wedge-shaped region be $r$ and the angle be $\theta$, then the stress field is expressed as:

Here $K$ is a coefficient related to stress intensity, $f_{ij}(\theta)$ is an angle-dependent function, and $\lambda$ is the Williams eigenvalue. This $\lambda$ is determined depending on the wedge opening angle $2\alpha$. The key point is that when $\lambda < 1$, $\sigma \to \infty$ as $r \to 0$.

Solving the Williams eigenvalue equation gives $\lambda \approx 0.545$ for an internal angle of 270° (reentrant corner). That means stress diverges on the order of $\sigma \propto r^{-0.455}$. For a crack tip (internal angle 360°), it's the familiar $\lambda = 0.5$ and $\sigma \propto r^{-1/2}$.

Exactly. So in practice, we avoid singularities by adding a fillet R (rounding) to eliminate the corner. However, even if the CAD model has a fillet, if the mesh is too coarse to adequately represent the corner, a numerical singularity may remain, so caution is needed.

To be precise, stresses at positions away from the singularity do converge properly. The problem is the maximum stress right at the singularity. Let's look at the standard FEM error estimate:

Here $k$ is the element polynomial order, and $\lambda$ is the Williams eigenvalue. For a smooth solution, $\lambda \geq k$ and we get a convergence rate of $h^k$. However, with a singularity, $\lambda < 1 < k$, so the convergence rate is limited to $h^\lambda$. Even using quadratic elements, it only reaches about $h^{0.545}$, not $h^2$.

It's not meaningless for the global energy norm error, just that improvement is slow. However, the stress value at the singularity itself does not converge even if you increase the element order or refine the mesh—this is a fundamental limitation. It's not a bug in FEM, but an inherent property of the physical model (a sharp corner in a linear elastic body).

Exactly. So "which stress to evaluate" is the essence of the singularity problem. Looking at a maximum stress contour plot and reflexively judging "The red area is dangerous!" is completely wrong for models with singularities.

Let me list some typical ones:

• Reentrant Corner — L-shaped bracket, T-joint corner without fillet. The most frequently encountered in practice.
• Crack Tip — Intentionally introduced in fracture mechanics models. $\lambda = 0.5$, the strongest singularity.
• Edge of Dissimilar Material Interface — Junction edge of materials with different Young's modulus. Called a bi-material singularity.
• Point of Concentrated Load Application — Point load/line load theoretically has zero contact area → infinite stress.
• Corner of Constraint Condition — Corner at the boundary between a fixed end and a free end. Creates a numerical singularity.
• V-Notch — Weld bead termination, keyway bottom, etc.

Yes. If you replace a bolt seating surface or bearing contact surface with a "concentrated load on a single node", the stress near that node is physically meaningless. At a minimum, you must distribute the load over an area or exclude stresses near the load application point from evaluation. That's the golden rule.

Good question. In fracture mechanics, we evaluate using "stress intensity factor $K$" or "$J$-integral", not the "stress value" itself. However, as a technique to accurately obtain $K$ with FEM, there are quarter-point elements (1/4-point elements).

A method devised by Barsoum (1976). For an 8-node quadrilateral quadratic element, shift the mid-side node on the crack tip side to the quarter-point position along the edge. Then a $\sqrt{r}$ term naturally appears in the element displacement field, allowing accurate capture of the $r^{-1/2}$ singularity of the crack tip at the element level.

Abaqus and Ansys have built-in functionality to place these quarter-point elements around crack tips. Combined with the J-integral or contour integral methods mentioned later, this allows high-precision determination of $K_I, K_{II}, K_{III}$.

hp-refinement is effective. Refine the mesh geometrically (geometric grading) towards the singularity while simultaneously increasing the element polynomial order. According to research by Babuška et al., this combination can achieve exponential convergence:

Here $N$ is the number of degrees of freedom. While standard h-refinement only converges algebraically ($N^{-\lambda/d}$), hp-refinement converges exponentially, making it vastly more efficient. However, implementation is complex and only limited solvers support it. StressCheck and p-FEM based solvers are representative examples.

There are mainly three methods:

1. Displacement Extrapolation Method — Extrapolate using the relation $K = \lim_{r \to 0} \sigma\sqrt{2\pi r}$ from nodal displacements near the crack tip. Simple but highly mesh-dependent.
2. J-Integral Method (Equivalent Domain Integral) — Calculate the energy release rate $G$ along a path (contour) surrounding the crack tip. $K_I = \sqrt{E' G}$ ($E'$ differs for plane stress/strain). Path-independent and high accuracy.
3. Interaction Integral Method — Separately extract $K_I, K_{II}, K_{III}$ using the interaction between the J-integral of an auxiliary field and the real field. Essential for mixed-mode problems.

Exactly. So we don't use the stress value right at the singularity, but find $J$ using the average of multiple slightly distant contours. In Abaqus, *CONTOUR INTEGRAL automatically calculates 5-10 contours. Outer contours are less affected by mesh and more stable.

The industry standard is stress linearization as defined by ASME Boiler & Pressure Vessel Code Section VIII Division 2. This is a method to decompose the stress distribution through the thickness into "membrane stress $\sigma_m$" and "bending stress $\sigma_b$".

Here $t$ is the plate thickness, $x$ is the coordinate through the thickness. Membrane stress is the component pulling the cross-section uniformly, bending stress is the component with opposite signs on the front and back surfaces. The peak stress $\sigma_p$ caused by the singularity (the residual nonlinear distribution after linearization) is excluded from primary evaluation.

It is ignored for "primary stress evaluation". However, peak stress becomes important for fatigue evaluation. ASME Div.2 Part 5.5 uses full-range stress (membrane+bending+peak) for fatigue analysis. At this point, the peak stress at the singularity is physically meaningless, so the correct approach is to switch to a sub-model with a fillet R and evaluate the actual stress concentration.