Linear Stress Recovery Methods in Elasticity

• Linear stress recovery methods are post-processing techniques that produce statically admissible and superconvergent stress approximations from raw finite element solutions.
• They employ patch-based least-squares and polynomial-preserving recovery operators to achieve second-order convergence (O(h²)) even when the underlying solution converges only at first order.
• These methods are practical for adaptive mesh refinement and error estimation, validated across regular and distorted meshes in hybrid-stress finite element frameworks.

Linear stress recovery methods are a foundational class of post-processing procedures in computational elasticity, devised to produce statically admissible and superconvergent approximations to physical stress fields from a computed displacement or stress solution. These methods are central to a posteriori error estimation, adaptive mesh refinement, and accurate physical interpretation of computed results, especially where raw numerical stress fields exhibit suboptimal convergence or discontinuities. Theoretical developments and practical algorithms span primal finite elements, hybrid-stress elements, smoothed FEM, and domain decomposition, with significant impact on both mathematical understanding and engineering reliability.

1. Formulation in the Hybrid-Stress Finite Element Method

A paradigmatic instance of linear stress recovery is the methodology introduced for the 4-node hybrid-stress quadrilateral (Pian–Sumihara) element. The canonical setting seeks $(\sigma,u) \in E \times V$, where $$E = \{\tau \in L^2(\Omega; \mathbb{R}^{2 \times 2}_{\text{sym}}) : \operatorname{tr} \tau = 0\}$$ and $V = \{v \in (H_0^1(\Omega))^2\}$. The hybrid formulation solves

$$a(\sigma,\tau) + b(\tau,u) = 0, \quad b(\sigma,v) = F(v),$$

with $a$ and $b$ standard elasticity bilinear forms. The Pian–Sumihara stress space $E_h$ on each element $K$ is constructed using a 5-parameter mode: $$\hat\sigma_h = \begin{pmatrix} 1 & 0 & 0 & \eta & \xi \ 0 & 1 & 0 & \frac{b_2}{b_1}\xi & \frac{a_1}{a_2}\eta \ 0 & 0 & 1 & \xi & \eta \end{pmatrix} \beta,$$ where $\beta \in \mathbb{R}^5$, with local coordinates $(\xi,\eta)$. The global stress solution is piecewise but generally discontinuous between elements (Bai et al., 2015).

2. Vertex Patch-Based Least-Squares Recovery Operator

The core stress recovery operator $R_h$ maps the raw discontinuous fem stress $\sigma_h$ to a globally continuous, piecewise-bi-linear field $M_h$. For each mesh vertex $Z_i$ with adjacent patch $\omega_i$, one solves, for each component,

$$p_i(x,y) = \alpha^i_1 + \alpha^i_2 x + \alpha^i_3 y,$$

minimizing

$$J(p) = \sum_{K \subset \omega_i} \int_K (p(x,y) - \sigma_h^\bullet(x,y))^2 dx\,dy.$$

The normal matrix and right-hand side are assembled over the patch, and $R_h\sigma_h^\bullet (Z_i)$ is set to $p_i(Z_i)$. Extension to the full domain uses the standard bi-linear shape functions that define $M_h$. This construction ensures $R_h$ is linear, bounded, and polynomial-preserving on each patch (Bai et al., 2015).

3. Mesh Distortion, Superconvergence, and Error Estimation

Analytical superconvergence for the recovered stress depends on a mesh distortion criterion: for each quadrilateral $K$, the distance $d_K$ between the diagonal midpoints satisfies $d_K = O(h_K^\alpha)$, with $\alpha>0$. Under this diagonal condition and mild “neighboring” regularity, superconvergence of the recovered stress is demonstrated: $$\|\sigma - R_h\sigma_h\|_{L^2(\Omega)} \leq C \left( h^{1+\alpha} \|\sigma\|_{H^1} + h^2 \|\sigma\|_{H^2} \right) = O(h^{1 + \min\{\alpha,1\}}).$$ Numerical experiments confirm that under practical mesh sequences with $\alpha=1$ (e.g., newest-vertex bisection), full second-order convergence $\|\sigma - R_h\sigma_h\| = O(h^2)$ is observed, despite the underlying solution converging only as $O(h)$ (Bai et al., 2015).

For a posteriori error control, define

$$\eta_\sigma = \| R_h\sigma_h - \sigma_h \|_{L^2(\Omega)}.$$

Then

$$|\| \sigma - \sigma_h \| - \eta_\sigma| \leq \| \sigma - R_h\sigma_h \|,$$

so $$\eta_\sigma = \| \sigma - \sigma_h \| [1 + O(h^{1 + \min\{\alpha,1\}})]$$, i.e., the estimator is asymptotically exact.

4. Recovery of Displacement Gradients and Dual Superconvergence

Analogous polynomial-preserving recovery (PPR) is employed for gradients of the displacement field. At each node, fit a quadratic polynomial in a least-squares sense to neighboring displacement nodal values on the patch and set the recovered gradient as its derivative at the node, extended globally by bilinear interpolation. This achieves the same superconvergent rate under identical mesh conditions: $$\| \nabla u - G_hu_h \|_{L^2(\Omega)} \leq C \left(h^{1+\alpha} \|u\|_{H^3} + h^2\|\sigma\|_{H^2}\right).$$ The error estimator

$\eta_u = \|G_hu_h - \nabla u_h\|_{L^2(\Omega)}$

is likewise asymptotically exact, with estimator-effectivity tending to unity as $h \to 0$ (Bai et al., 2015).

5. Numerical Investigation and Practical Impact

Extensive tests on both regular and highly irregular quadrilateral meshes (unit square, rectangles with large aspect ratio, and near-incompressible Poisson ratio up to 0.4999) confirm that the stress and gradient recoveries achieve true second-order convergence even where the base FE solution is merely first-order:

• $\|\sigma_h-\sigma\| = O(h)$,
• $\|\sigma - R_h \sigma_h\| = O(h^2)$,
• $\eta_\sigma / \|\sigma - \sigma_h\| \to 1$ as $h \to 0$.

This robust superconvergence is attained using only local patch-based linear least-squares fits, without requiring data beyond neighboring element stresses. The methodology is thus practical for implementation in standard finite element codes and supports reliable adaptive mesh refinement cycles (Bai et al., 2015).

6. Context within Broader Stress Recovery Methodologies

The linear least-squares patch recovery operator employed here is a specific instance within a constellation of recovery-based approaches, which also include Zienkiewicz–Zhu Simple Averaging, Superconvergent Patch Recovery (SPR), and Polynomial-Preserving Recovery (PPR) methods (Guo et al., 2024). All these approaches fundamentally aim to construct a continuous, higher-regularity approximation to the stress field converging at a rate faster than that of the underlying FE solution, under mesh regularity conditions. Recovery-based estimators are widely used for reliable a posteriori error control and mesh adaptation in computational mechanics.

7. Limitations and Prospects

The theoretical superconvergence relies critically on the mesh distortion parameter $\alpha$ and proper patch construction around vertices. For highly distorted, non-quadrilateral meshes, convergence rates may be suboptimal unless mesh regularity is enforced. The method is tailored to hybrid-stress elements; for higher-order elements or other mixed methods, alternative local recovery constructions or equilibrium-constrained recovery techniques may be necessary.

Linear stress recovery methods, underscored by the patchwise least-squares strategy of Bai, Wu, and Xie, are a mature and robust toolset for post-processing, error estimation, and adaptive computation in linear elasticity (Bai et al., 2015).