Self-Stabilized Virtual Elements Overcome Artificial Hardening in 2D Plasticity Computations

Computational engineers Ossama Firari, Massimiliano Cremonesi, and Umberto Perego have formulated a variationally consistent numerical framework based on the Virtual Element Method to simulate two-dimensional elastoplastic deformation without suffering from spurious mesh locking or unphysical artificial hardening. The research, published in Advanced Modeling and Simulation in Engineering Sciences, establishes a mathematical formulation that successfully tracks plastic failure limits across highly distorted polygonal meshes, including both convex and non-convex quadrilaterals and pentagons. By extending the classical mixed Hu-Washizu variational functional of elasticity to encompass the energy contributions governing irreversible finite-step elastoplastic transformations, the authors demonstrated that all fundamental governing relations emerge directly from the functional’s stationarity conditions. These deduced relations encompass static equilibrium, elastic constitutive state equations, geometric strain-displacement compatibility, internal variable evolution laws, and the complementary Kuhn-Tucker loading-unloading conditions.

The Virtual Element Method represents a modern generalization of the conventional Finite Element Method, providing structural analysts with the geometric freedom to discretize complex two-dimensional domains into arbitrary polygons with any number of edges and arbitrary convex or non-convex profiles. Despite this topological versatility, applying the method to nonlinear elastoplasticity has historically encountered severe computational hurdles related to numerical stability and volumetric locking. Standard formulations typically require adding artificial, parameter-dependent stabilization stiffness matrices to suppress zero-energy hourglass deformation modes. Furthermore, conventional low-order elements frequently exhibit pathological volumetric locking when plastic deformation approaches the theoretical collapse threshold under unconstrained plastic flow conditions, manifesting as an incorrect, artificially steep hardening response rather than the expected constant-load limit state.

To overcome these limitations without tuning artificial stabilization parameters, the investigators established a self-stabilized formulation that matches displacement interpolations along element boundaries with independent internal strain approximations of identical polynomial order. For first-order virtual elements where displacement fields vary linearly along element perimeters, the authors adopted linear, non-divergence-free strain representations characterized by either seven or nine independent parameters. Because the resulting stiffness matrix satisfies the dimensional stability requirement where the difference between displacement degrees of freedom and strain parameters remains at or below three, hourglass kinematic modes are suppressed entirely from the internal energy construction, providing rigorous stabilization directly within the variational architecture.

A critical analytical development in the formulation involves resolving domain integrals that govern the element compatibility matrix, which naturally contain internal displacement fields whose explicit algebraic expressions are intentionally undefined within virtual elements. By executing an integration by parts over the element domain, the compatibility integration separates into a fully computable boundary integral governed by known boundary polynomials and a remaining internal domain integral. The authors showed that evaluating this remaining internal domain term requires incorporating the first two internal moments of the unknown displacement shape functions as non-nodal degrees of freedom. For a four-node quadrilateral element, incorporating these two moment degrees of freedom increases the displacement unknown count from eight to ten, whereas a five-node pentagon expands from ten to twelve displacement unknowns.

The incremental elastoplastic boundary value problem is solved computationally using a backward-difference return-mapping scheme integrated with a global Newton-Raphson iterative algorithm. In this numerical workflow, trial elastic stress states are systematically evaluated at internal integration points defined across triangular sub-elements. If the yield criterion is exceeded, local stress tensors, plastic strain increments, and kinematic internal variables are iteratively corrected, simultaneously yielding the consistent elastoplastic tangent matrix. The authors mathematically derived the exact closed-form representation of this tangent matrix for the mixed virtual element discretization, ensuring an asymptotic second-order quadratic convergence rate during global equilibrium iterations until residual force norms satisfy prescribed convergence tolerances.

The performance of the self-stabilized elements was evaluated across three classic benchmark problems involving perfectly plastic materials known to cause severe pathological locking in standard computational formulations. The first benchmark evaluated the ultimate bearing capacity of a rigid, frictionless strip footing resting on a weightless soil medium governed by the Drucker-Prager yield criterion. The domain, discretized into symmetric half-space models, utilized an elastic Young’s modulus of 207 megapascals, a Poisson’s ratio of 0.3, a cohesion of 0.069 megapascals, and an internal friction angle of 20 degrees. The theoretical limit pressure established by Prandtl’s limit equilibrium theory yields an exact collapse value of 1.024 megapascals. The simulations subjected the foundation to an incremental downward displacement of 0.065 meters across meshes featuring 325 convex quadrilaterals, 368 mixed convex and non-convex quadrilaterals, and 368 mixed pentagonal elements.

When analyzed with standard stabilized virtual elements or conventional fully integrated four-node isoparametric finite elements, the strip footing simulations exhibited severe volumetric locking characterized by continuous unphysical load escalation. Similarly, self-stabilized quadrilateral formulations omitting the internal moment degrees of freedom failed to track the limit state. In contrast, the self-stabilized virtual elements incorporating the two internal moment degrees of freedom successfully captured the theoretical plastic collapse limit without exhibiting artificial hardening. These enhanced elements achieved close agreement with benchmark reference solutions computed using ABAQUS CPE4RH elements, which incorporate reduced integration, hourglass control, and constant independent pressure assumptions, thereby demonstrating the decisive role of displacement moments in eliminating locking behavior across irregular grids.

The second numerical evaluation investigated a plane strain double-edge-notched specimen subjected to tensile deformation governed by von Mises perfect plasticity. The geometry, defined by length-to-width and width-to-ligament ratios of 3 and 10, was tested under prescribed vertical displacements up to 0.3 millimeters at the top boundary. Material properties included an elastic modulus of 70 megapascals, a Poisson’s ratio of 0.3, and a yield stress of 0.243 megapascals, corresponding to a theoretical Prandtl kinematic limit net stress of approximately 2.97 times the yield strength. Testing regular quadrilateral meshes, non-convex quadrilateral meshes, and mixed pentagonal arrangements within the plastic ligament observation zone confirmed that standard virtual elements developed steep linear hardening following initial yield. Conversely, both the seven-parameter and nine-parameter self-stabilized formulations with internal moments tracked the theoretical net limit stress plateau precisely while accurately mapping plastic strain concentrations localized around the notch tip.

The third benchmark examined an annular circular ring with an inner radius of 10 millimeters and a thickness of 1 millimeter subjected to diametrically opposed vertical displacements under displacement control. Operating under predominantly flexural plastic action localized along narrow bands near the axes of symmetry, the structure modeled a von Mises material having a Young’s modulus of 1,000,000 megapascals, a Poisson’s ratio of 0.25, and a yield stress of 40 megapascals, associated with an exact kinematic limit load of 4.404 newtons. Numerical evaluations conducted across regular meshes, non-convex quadrilateral meshes, and pentagonal arrangements revealed that while self-stabilized formulations lacking internal moments consistently locked across all mesh topologies, the moment-enriched virtual elements converged accurately to the 4.404 newton plastic collapse plateau without numerical anomalies.

The authors noted specific theoretical boundaries regarding stability proofs and structural scope across their findings. Although the mathematical stability of the consistent stiffness matrix has been proven for linear quadrilateral elements, equivalent mathematical stability proofs for pentagons and polygons with higher edge counts currently rely on necessary dimensional inequalities validated numerically rather than exhaustive analytical theorems. The published investigation remains focused on small-displacement kinematics in two-dimensional plane strain settings. Nonetheless, the findings establish that integrating internal displacement moment degrees of freedom into mixed variational virtual element frameworks provides an accurate, locking-free methodology capable of handling extreme element distortions and arbitrary polygonal geometries in nonlinear solid mechanics.

Subject of Research: Self-stabilized virtual element methods for two-dimensional locking-free elastoplastic computational modeling

Article Title: Variationally consistent self-stabilized Virtual Elements for 2D locking-free elastoplasticity

Article References: Firari, O., Cremonesi, M., & Perego, U. (2025). Variationally consistent self-stabilized Virtual Elements for 2D locking-free elastoplasticity. Advanced Modeling and Simulation in Engineering Sciences, 12(1), Article 13. https://doi.org/10.1186/s40323-025-00297-3

Image Credits: AI Generated

DOI: 10.1186/s40323-025-00297-3

Keywords: Virtual Element Method, elastoplasticity, computational mechanics, finite element method, hourglass stabilization, volumetric locking, polygonal meshes, Hu-Washizu functional, strip footing, Drucker-Prager plasticity, von Mises plasticity, internal moment degrees of freedom, element distortion, plastic collapse load, return-mapping algorithm, Newton-Raphson iteration, non-convex elements