Nitsche methods for constrained problems in mechanics

Abstract: We present guidelines for deriving new Nitsche Finite Element Methods to enforce equality and inequality constraints that act on the value of the unknown mechanical quality. We first formulate the problem as a stabilized finite element method for the saddle point formulation where a Lagrange multiplier enforces the underlying constraint. The Nitsche method is then presented in a general minimization form, suitable for nonlinear finite element methods and allowing straightforward computational implementation with automatic differentation. To validate these ideas, we present Nitsche formulations for a range of problems in solid mechanics and give numerical evidence of the convergence rates of the Nitsche method.

• The paper develops a stabilized saddle-point framework that eliminates Lagrange multipliers elementwise and yields a general minimization template for equality and inequality constraints.
• The framework derives methods for membrane, solid, and plate contact plus Kirchhoff plate boundary inequalities, with experiments showing expected optimal convergence rates.
• The approach can improve conditioning over penalty methods, but rigorous stability and error analyses for the new formulations remain open research problems.

Overview

This paper presents a systematic framework for deriving Nitsche-type finite element methods for mechanical problems subject to equality and inequality constraints acting on the values of the unknown field. Rather than treating the Nitsche method as a consistency correction to the penalty method—the prevailing modern viewpoint—the authors build on Stenberg's reinterpretation of Nitsche's method as an element-wise elimination of the Lagrange multiplier from a residual-stabilized saddle point formulation. From this starting point, they distill a general minimization form applicable to arbitrary constrained minimization problems, and use it to derive and implement several novel methods: two-membrane contact, membrane–solid contact, plate–plate contact, and a Kirchhoff plate with an inequality boundary condition. Numerical experiments confirm optimal convergence rates in each case.

The stabilized saddle point viewpoint

The paper begins with the classical model problem $-\nabla\cdot(\kappa\nabla u)=f$ with Dirichlet condition $u=g$, discretized in $V_h\subset H^1(\Omega)$. Following Stenberg [stenberg1995], the Nitsche method is obtained by stabilizing the mixed formulation with Lagrange multiplier $\lambda=\kappa\,\partial u/\partial n$ via Barbosa–Hughes residual stabilization:

$$\mathcal{L}_h(v_h,\mu_h) = \mathcal{L}(v_h,\mu_h) - \int_{\partial\Omega}\frac{\gamma}{2}\Big(\mu_h - \kappa\frac{\partial v_h}{\partial n}\Big)^2\,\mathrm{d}s,$$

with $\gamma|_K = \alpha h_K/\kappa$. Stability holds for any conforming pair $(V_h,Q_h)$ provided $0<\alpha<C_I$, where $C_I$ is an inverse-estimate constant—an important robustness property compared to inf–sup-constrained mixed methods. Element-wise elimination of $\lambda_h$ (exact when $u=g$0 contains discontinuous polynomials of degree at least that of $u=g$1, so that the $u=g$2 projection is the identity) recovers precisely the classical symmetric Nitsche functional

$u=g$3

The authors contrast this with the penalty method, which effectively enforces a Robin condition $u=g$4; driving $u=g$5 yields ill-conditioned systems, whereas the Nitsche form retains consistency without this deterioration.

Inequality constraints: Signorini's problem

For the scalar Signorini problem ($u=g$6 on $u=g$7), the discrete multiplier is the orthogonal $u=g$8 projection of the residual expression onto nonnegative functions. Choosing $u=g$9 makes this projection explicit via the maximum operator,

$V_h\subset H^1(\Omega)$0

and substitution into the stabilized Lagrangian produces the inequality Nitsche functional, identical to the equality case except for the positive-part operator $V_h\subset H^1(\Omega)$1. This observation is the key structural insight: equality constraints are recovered simply by dropping the maximum operator, and the penalty method is recovered by formally setting $V_h\subset H^1(\Omega)$2.

A general template

The general form of the method is stated as: find $V_h\subset H^1(\Omega)$3 minimizing

$V_h\subset H^1(\Omega)$4

where three ingredients must be supplied for each new problem:

1. Constraint: $V_h\subset H^1(\Omega)$5 on the candidate active set $V_h\subset H^1(\Omega)$6.
2. Continuous Lagrange multiplier: $V_h\subset H^1(\Omega)$7, obtained from the strong form of the continuous saddle point problem.
3. Scaling: $V_h\subset H^1(\Omega)$8, chosen so that $V_h\subset H^1(\Omega)$9 with mesh- and material-independent proportionality—physically, $\lambda=\kappa\,\partial u/\partial n$0 converts displacement units into force units, and mathematically it is required for uniform stability via inverse estimates.

For two-body contact, prior error analysis indicates that $\lambda=\kappa\,\partial u/\partial n$1 and hence $\lambda=\kappa\,\partial u/\partial n$2 should be defined on the less stiff side of the interface. The functional is posed directly as a minimization problem, which suits nonlinear problems: Newton linearization is performed on the energy itself, with first and second Gateaux derivatives computed automatically via JAX within the scikit-fem assembly framework.

Existing methods recovered

Two known formulations fit the template. For the membrane obstacle problem ($\lambda=\kappa\,\partial u/\partial n$3 in $\lambda=\kappa\,\partial u/\partial n$4), the multiplier is $\lambda=\kappa\,\partial u/\partial n$5 with scaling $\lambda=\kappa\,\partial u/\partial n$6. For two-body linearized elastic contact, the constraint is $\lambda=\kappa\,\partial u/\partial n$7 on the contact interface, the multiplier is the normal traction $\lambda=\kappa\,\partial u/\partial n$8, and $\lambda=\kappa\,\partial u/\partial n$9; the resulting functional coincides with the symmetric variant of Chouly, Hild, and Renard [chouly2015symmetric].

Novel methods and numerical validation

Two-membrane contact. Two Poisson membranes separated by gap $$\mathcal{L}_h(v_h,\mu_h) = \mathcal{L}(v_h,\mu_h) - \int_{\partial\Omega}\frac{\gamma}{2}\Big(\mu_h - \kappa\frac{\partial v_h}{\partial n}\Big)^2\,\mathrm{d}s,$$0 satisfy $$\mathcal{L}_h(v_h,\mu_h) = \mathcal{L}(v_h,\mu_h) - \int_{\partial\Omega}\frac{\gamma}{2}\Big(\mu_h - \kappa\frac{\partial v_h}{\partial n}\Big)^2\,\mathrm{d}s,$$1 in $$\mathcal{L}_h(v_h,\mu_h) = \mathcal{L}(v_h,\mu_h) - \int_{\partial\Omega}\frac{\gamma}{2}\Big(\mu_h - \kappa\frac{\partial v_h}{\partial n}\Big)^2\,\mathrm{d}s,$$2, with contact pressure $$\mathcal{L}_h(v_h,\mu_h) = \mathcal{L}(v_h,\mu_h) - \int_{\partial\Omega}\frac{\gamma}{2}\Big(\mu_h - \kappa\frac{\partial v_h}{\partial n}\Big)^2\,\mathrm{d}s,$$3 derived from the less stiff membrane. Because $$\mathcal{L}_h(v_h,\mu_h) = \mathcal{L}(v_h,\mu_h) - \int_{\partial\Omega}\frac{\gamma}{2}\Big(\mu_h - \kappa\frac{\partial v_h}{\partial n}\Big)^2\,\mathrm{d}s,$$4 appears undifferentiated while $$\mathcal{L}_h(v_h,\mu_h) = \mathcal{L}(v_h,\mu_h) - \int_{\partial\Omega}\frac{\gamma}{2}\Big(\mu_h - \kappa\frac{\partial v_h}{\partial n}\Big)^2\,\mathrm{d}s,$$5 involves second derivatives, dimensional consistency forces $$\mathcal{L}_h(v_h,\mu_h) = \mathcal{L}(v_h,\mu_h) - \int_{\partial\Omega}\frac{\gamma}{2}\Big(\mu_h - \kappa\frac{\partial v_h}{\partial n}\Big)^2\,\mathrm{d}s,$$6. With linear elements (for which the element-wise Laplacian vanishes), the observed convergence rate is linear in the $$\mathcal{L}_h(v_h,\mu_h) = \mathcal{L}(v_h,\mu_h) - \int_{\partial\Omega}\frac{\gamma}{2}\Big(\mu_h - \kappa\frac{\partial v_h}{\partial n}\Big)^2\,\mathrm{d}s,$$7 norm, matching theory. A notable comparative result concerns conditioning: with quadratic elements, the penalty method requires $$\mathcal{L}_h(v_h,\mu_h) = \mathcal{L}(v_h,\mu_h) - \int_{\partial\Omega}\frac{\gamma}{2}\Big(\mu_h - \kappa\frac{\partial v_h}{\partial n}\Big)^2\,\mathrm{d}s,$$8 to retain optimal convergence, and consequently exhibits consistently larger Jacobian condition numbers and slower Newton convergence than the Nitsche variant with $$\mathcal{L}_h(v_h,\mu_h) = \mathcal{L}(v_h,\mu_h) - \int_{\partial\Omega}\frac{\gamma}{2}\Big(\mu_h - \kappa\frac{\partial v_h}{\partial n}\Big)^2\,\mathrm{d}s,$$9.

Membrane against elastic solid. A Poisson membrane contacts a linear elastic cube across a gap; here $\gamma|_K = \alpha h_K/\kappa$0 and $\gamma|_K = \alpha h_K/\kappa$1. Linear hexahedral/quadrilateral elements yield the expected linear $\gamma|_K = \alpha h_K/\kappa$2 convergence rate.

Plate against plate. With biharmonic energies, the multiplier is $\gamma|_K = \alpha h_K/\kappa$3 and $\gamma|_K = \alpha h_K/\kappa$4. Using nonconforming Morley elements (again rendering the element-wise biharmonic term zero), the method attains the theoretical linear rate in the $\gamma|_K = \alpha h_K/\kappa$5 norm.

Kirchhoff plate with inequality boundary condition. For the plate corners problem with constraint $\gamma|_K = \alpha h_K/\kappa$6 on $\gamma|_K = \alpha h_K/\kappa$7, the multiplier is the Kirchhoff shear force $\gamma|_K = \alpha h_K/\kappa$8 and $\gamma|_K = \alpha h_K/\kappa$9. With Bogner–Fox–Schmit elements, quadratic convergence is observed, consistent with the element's polynomial degree.

Across all four problems, convergence rates were estimated from differences of solutions on successively refined meshes ($(V_h,Q_h)$0 vs. $(V_h,Q_h)$1), which bounds the true error rate under the assumption of asymptotic behavior $(V_h,Q_h)$2.

Limitations and open questions

The paper is explicitly a methodology and numerics contribution: proofs of stability and a priori error estimates for the newly proposed formulations are deferred to future work, so the optimality of the observed rates rests on numerical evidence rather than analysis. Several choices also remain heuristic or empirically grounded rather than proven: the rule that stabilization should be placed on the less stiff body (or the body with smaller elements) is adopted from prior error analyses; the stabilization constant $(V_h,Q_h)$3 is fixed by hand in each experiment (e.g., $(V_h,Q_h)$4 for the contact problems, $(V_h,Q_h)$5 for the plate corner problem); and the convergence-rate estimation via successive-mesh differences assumes the asymptotic regime has been reached. Whether the general functional admits uniform stability and optimal error estimates for arbitrary constraints—and how sensitive the results are to $(V_h,Q_h)$6—remain open questions.

Conclusion

The paper reframes the derivation of Nitsche methods for constrained mechanics around a stabilized saddle point structure, yielding a compact minimization template parameterized by the constraint $(V_h,Q_h)$7, the continuous multiplier $(V_h,Q_h)$8, and the scaling $(V_h,Q_h)$9. The template reproduces existing methods for obstacle and elastic contact problems and generates new ones for two-membrane, membrane–solid, and plate–plate contact, plus an inequality boundary condition for Kirchhoff plates, all implemented through automatic differentiation and validated at their theoretical convergence rates. Rigorous stability and error analysis of the new formulations constitutes the principal outstanding task.