Adaptive finite elements for obstacle problems

We summarise three applications of the obstacle problem to membrane contact, elastoplastic torsion and cavitation modelling, and show how the resulting models can be solved using mixed finite elements. It is challenging to construct fixed computational meshes for any inequality-constrained problem because the coincidence set has an unknown shape. Consequently, we demonstrate how $h$-adaptivity can be used to resolve the unknown coincidence set. We demonstrate some practical challenges that must be overcome in the application of the adaptive method.

💡 Research Summary

This paper investigates the numerical solution of obstacle problems that arise in three classical engineering applications: membrane‑contact, elastoplastic torsion, and hydrodynamic lubrication with cavitation. The authors begin by formulating the classical obstacle problem on a two‑dimensional domain Ω as a variational inequality with the Karush‑Kuhn‑Tucker conditions
‑Δu ≥ f, u ≥ g, (‑Δu − f)(u − g) = 0, u|∂Ω = 0.
The unknown coincidence set Ω_C where u = g (contact) is the main difficulty because its shape depends on the solution itself, making a fixed mesh inefficient. To overcome this, the paper adopts an h‑adaptive strategy combined with a mixed finite‑element discretisation: a bubble‑enriched quadratic space P₂ ⊕ B₃ for the primal variable u and a piecewise‑constant space P₀ for the Lagrange multiplier λ (interpreted as contact pressure or a cavitation indicator).

Each application is recast as an obstacle problem. In the membrane case, the membrane tension κ leads to –∇·(κ∇u) ≥ f with the unilateral constraint u ≥ g (rigid obstacle). In the torsion problem, the stress function φ satisfies a gradient bound |∇φ| ≤ τ√3; after a transformation ψ = –φ the problem becomes a standard obstacle formulation with a distance‑to‑boundary obstacle. In the lubrication case, the Reynolds equation is coupled with the Swift–Stieber cavitation condition p ≥ p_cav, turning the pressure field into an obstacle problem where the cavitation region is Ω_C.

The continuous variational formulation introduces λ as a Lagrange multiplier, yielding a saddle‑point problem: find (u,λ)∈V×Λ such that
(∇u,∇v) − ⟨λ, v⟩ = (f, v) ∀v∈V,
⟨μ − λ, u − g⟩ ≥ 0 ∀μ∈Λ.
Here V = H₀¹(Ω) and Λ is the cone of non‑negative functionals. Discretisation replaces λ∈H^{‑1}(Ω) by λ_T∈Q_T (piecewise constants) and restricts the positivity cone to Λ_T. The mixed spaces (V_T, Q_T) satisfy a discrete inf‑sup condition, guaranteeing stability.

A key contribution is an a‑posteriori error estimator (Theorem 1) that bounds the combined H¹‑error of u and H^{‑1}‑error of λ by elementwise residuals: interior residual η_K, edge jump η_∂K, and a contact‑constraint residual η_{C,K}. This estimator drives the adaptive mesh refinement: elements with large η_K are marked and refined, concentrating degrees of freedom near the unknown free boundary.

Because the inequality constraints make the problem nonlinear, the authors employ a primal‑dual active‑set (PDAS) method, which is mathematically equivalent to a semismooth Newton iteration. The method projects the discrete multiplier onto the non‑negative cone via λ_T = (λ_T − π_T(u_T − g))⁺, where π_T is the L²‑projection onto Q_T. Consequently, each element either has λ_T = 0 (inactive) or satisfies π_T(u_T − g) = 0 (active). The algorithm iterates: (1) identify the active set A where λ_T − π_T(u_T − g) > 0, (2) enforce λ_T = 0 on the complement, (3) solve the resulting linear mixed system for (u_T,λ_T) on the reduced space, and (4) check convergence via the relative change of λ_T. The authors provide a compact Python implementation using the scikit‑fem library, illustrating mesh creation, assembly of bilinear and linear forms, active‑set selection, and solution of the condensed linear system.

Implementation details follow the classic SOLVE‑ESTIMATE‑MARK‑REFINE loop. After each PDAS solve, the error estimator is evaluated, a Dörfler‑type marking selects a subset of elements, and uniform refinement (or newest‑vertex bisection) creates a new mesh. The process repeats until a prescribed tolerance is reached. Numerical experiments for all three applications demonstrate that the adaptive scheme accurately captures the free boundary, yields optimal (or even super‑convergent) rates in the energy norm, and significantly reduces the number of degrees of freedom compared with uniform refinement.

In conclusion, the paper presents a coherent framework that combines mixed finite‑element discretisation, rigorous a‑posteriori error estimation, and a semismooth Newton‑type active‑set algorithm to solve obstacle problems efficiently. The methodology is validated on realistic engineering models, showing that h‑adaptivity is essential for resolving the unknown coincidence set and that the bubble‑enriched P₂–P₀ pair provides both stability and accurate approximation of the Lagrange multiplier. This work offers a practical, theoretically sound toolbox for engineers and computational scientists dealing with inequality‑constrained PDEs.