Achieving Material Robustness via Symmetric Stress Finite Element Discretizations

Technical Summary: Achieving Material Robustness via Symmetric Stress Finite Element Discretizations

Problem Statement
In continuum mechanics, variational formulations such as the Hellinger–Reissner elasticity and velocity–stress formulations for incompressible flow involve symmetric stress tensors. The symmetry of the Cauchy stress tensor is physically tied to the conservation of angular momentum. When discretizing these problems using $H(\text{div})$-conforming finite elements, practitioners must choose between enforcing symmetry strongly (requiring discrete tensors to be pointwise symmetric) or weakly (enforcing symmetry via a Lagrange multiplier). While weak symmetry is often preferred in literature due to the perceived complexity of constructing strongly symmetric elements, this paper investigates the impact of this choice on numerical accuracy, particularly for materials with anisotropic constitutive laws.

Methodology
The authors develop a theoretical framework and a series of numerical benchmarks to compare strongly and weakly symmetric schemes.

1. Mathematical Formulation: The study considers a general linearized static problem for both fluids and solids, where the stress $\sigma$ relates to the strain $\varepsilon(u)$ and an external field term $\tilde{F}$ (representing anisotropy). The problem is cast as a saddle-point system.
2. Discretization Schemes:
   • Strongly Symmetric: The Johnson-Mercier-Křížek (JMK) scheme and the Hu-Zhang (HZk) scheme are used. These enforce $\sigma_h \in \Sigma^{\text{sym}}_h$ pointwise.
   • Weakly Symmetric: The PEERS and Arnold-Falk-Winther (AFWk) schemes are used. These introduce a Lagrange multiplier $\omega_h$ (representing the skew-symmetric part of the gradient) to enforce symmetry in a weak sense.
3. Benchmark Construction: The authors construct "manufactured solutions" where the exact stress is identically zero ($\sigma \equiv 0$) but the displacement and rotation fields are non-trivial. These solutions are derived from constitutive laws inspired by:
   • Linear isotropic solids.
   • Transversely isotropic solids (fiber-reinforced materials).
   • Polar fluids (liquid crystal polymer networks).
     The magnitude of the anisotropic forcing term is controlled by a scaling parameter $\delta$.
4. Theoretical Analysis: An abstract Hilbert space theory is developed for generic saddle-point problems. The authors define "structure-preserving" discretizations and derive a priori error estimates to explain the observed numerical behavior.

Key Contributions

• Definition of Material Robustness: The paper introduces the concept of "material robustness," defined as the property where a numerical scheme produces a zero-stress approximation for any exact zero-stress configuration, regardless of the constitutive law or the magnitude of the anisotropic forcing.
• Identification of Failure Modes: The authors demonstrate that weakly symmetric schemes can yield arbitrarily poor stress approximations (errors scaling with the forcing parameter $\delta$) even when the exact stress is zero. This occurs specifically when the exact rotation tensor $\omega = \text{skw}(\nabla u)$ does not belong to the discrete Lagrange multiplier space $\Xi_h$.
• Unifying Theory: A rigorous theory is presented showing that strongly symmetric schemes satisfying specific kernel inclusion properties ($\ker B_h \subseteq \ker B$) are inherently material robust. In contrast, weakly symmetric schemes lack this invariance unless the discrete rotation space contains the exact rotation tensor.
• Error Estimates: New a priori error estimates are derived for partially structure-preserving discretizations, explicitly showing how the stress error depends on the approximation quality of the rotation tensor in the weakly symmetric case.

Results
Numerical experiments on unstructured meshes confirm the theoretical predictions:

• Isotropic Case: For linear isotropic solids where the exact rotation is constant (and thus easily representable), both strong and weak schemes perform similarly, with errors dominated by round-off.
• Anisotropic Cases: For transversely isotropic solids and polar fluids, where the rotation tensor is complex:
  • Strongly Symmetric Schemes (JMK, HZ3): Consistently produce stress errors at the level of machine precision (zero stress state), independent of the scaling parameter $\delta$.
  • Weakly Symmetric Schemes (PEERS, AFWk): Produce stress errors that scale linearly with $\delta$. Even for very high-order schemes (e.g., AFW3), if the exact rotation tensor is not contained in the discrete space $\Xi_h$, the stress error remains large.
• Transient Behavior: A transient polar fluid example shows that while both schemes agree initially, the non-robust weakly symmetric scheme fails to relax to the correct zero-stress steady state, whereas the robust strongly symmetric scheme does.

Significance and Claims
The paper claims that the choice between strong and weak symmetry enforcement is not merely a matter of implementation complexity but fundamentally affects the physical fidelity of the numerical solution. The authors argue that enforcing symmetry strongly is crucial for material robustness, ensuring that numerical schemes do not generate spurious stresses in stress-free configurations, particularly in applications involving rotational invariants or anisotropic materials.

The work advocates for the broader adoption of strongly symmetric finite elements (such as JMK and Hu-Zhang), noting that these elements are increasingly available in open-source software like Firedrake and NGSolve. The authors conclude that for problems sensitive to angular momentum conservation, strongly symmetric discretizations provide a necessary guarantee of accuracy that weakly symmetric schemes cannot offer without restrictive conditions on the discrete spaces.