The first problem in contact mechanics.
In 1882 Heinrich Hertz, at age 24, found the pressure between two pressed elastic bodies. It is still the reference for every contact solver. We pressed a rigid cylinder into an elastic half-space with the active-set solver of physicsbase. We got the Hertz pressure ellipse back. The half-width and the peak pressure are within a few tenths of a percent. There is no penalty parameter and no fudge factor.
Contact is where the finite element method stops being linear. Two bodies come together. You do not know the region where they touch, or the pressure on it, until you solve the problem. The region depends on the deformation. The deformation depends on the pressure. Thus it is a boundary condition that you must find. Hertz solved the simplest version in closed form. That solution is the reference. If a contact code cannot reproduce Hertz, you cannot trust it on a harder problem.
What Hertz proved
Press a long rigid cylinder of radius \(R\) into a flat elastic half-space (plane strain, frictionless). The contact is a strip of half-width \(a\). The pressure on it follows a semi-ellipse:
\[ p(x) = p_0\sqrt{1-\left(\tfrac{x}{a}\right)^2}, \qquad a=\sqrt{\frac{4PR}{\pi E^\*}}, \qquad p_0=\frac{2P}{\pi a} \]
Here \(P\) is the load per unit length. \(E^\*\) is the contact modulus. For a rigid indenter on an elastic solid, \(E^\* = E/(1-\nu^2)\). Two things make this a real test. The pressure is not uniform. It is zero at the edges and has a peak in the middle. Also the contact width \(a\) is an output. You do not prescribe it. A solver must find both.
Contact without a penalty parameter
Most contact codes push the two surfaces together with a stiff penalty spring. They tune its stiffness until the penetration is "small enough." physicsbase uses an active-set method. It treats non-penetration as an exact inequality constraint, not an approximate force. The idea is the standard Signorini method:
Each node on the surface gets a rigid ceiling. This is the profile of the indenter, lowered by the indentation depth \(\delta\). A node at position \(x\) may not rise above \(-\delta + x^2/2R\). The solver then repeats a set of simple rules. First, it solves with the touching nodes pinned to their ceilings. Then it releases any pinned node whose reaction is now in tension, because the surface would have to pull the indenter, and frictionless contact cannot do this. Then it activates any free node that has gone through its ceiling. It repeats until the touching set stops changing. At convergence, each active node is exactly on the indenter. Each inactive node is below it. No contact force pulls. These are the exact Kuhn-Tucker conditions of the contact problem. There is no penalty stiffness to tune and no penetration to accept.
The result
This is a graded plane-strain mesh. It is finer near the contact zone. A rigid parabolic indenter presses it. The recovered pressures are on the Hertz ellipse across the whole contact. The contact patch that the solver found matches the theory:
| Quantity | physicsbase | Hertz (closed form) | error |
|---|---|---|---|
| Contact half-width a | 0.6050 | 0.6038 | +0.2% |
| Peak pressure p₀ | 0.03328 | 0.03318 | +0.3% |
| Pressure profile p(x)/p₀ | follows √(1−(x/a)²) | <2% | |
The half-width is within 0.2%. The peak pressure is within 0.3%. Each interior contact node is within about 2% of the elliptic profile. The small deviation is at the edge of contact, as it should be. There the pressure gradient is infinite, and any finite mesh must round the corner. The right panel shows that the constraint is exact. The elastic surface (blue) is on top of the rigid parabola (orange) inside the patch. It separates cleanly outside it.
Calling it
To an agent, contact is an analysis mode with a list of stops. Give each surface node a barrier. This is a limit on a displacement component and the rigid side. Then physicsbase returns the deformation and the contact force at each node that touched:
POST /v1/solve
{ "analysis": "contact",
"nodes": […], "materials": {…}, "elements": [ /* quad4, plane_strain */ ],
"supports": [ /* symmetry + base */ ],
"contacts": [
{"node":i,"dof":"uy","side":"upper","limit": /* indenter profile at x_i */}, … ] }
← displacements + per-node contact_force + the active contact set
Honest footnotes
This is frictionless, small-strain, linear-elastic contact against a rigid obstacle. These are the Hertz assumptions, so the comparison is fair. We approximate the half-space by a finite block. It is large enough (many contact-widths in each direction) that its boundaries do not matter. The contact edge is at the local element size, so the half-width steps with the mesh. We report the pressure profile and peak, which converge cleanly. We do not rely on the single edge node. The next steps are deformable-on-deformable contact, friction (Coulomb stick-slip), and the sub-surface shear peak at \(0.48a\), which is the real reason that rolling contacts fail. The active-set core here extends to all of them. This shows that a general-purpose engine, with an exact-constraint contact solver and no tuned penalty, reproduces the 1882 result that every contact code must match.
References
- H. Hertz, "Über die Berührung fester elastischer Körper," Journal für die reine und angewandte Mathematik 92 (1882) 156–171.
- K. L. Johnson, Contact Mechanics, Cambridge University Press (1985), ch. 4.
- P. Wriggers, Computational Contact Mechanics, 2nd ed., Springer (2006).
- A. Signorini, "Sopra alcune questioni di elastostatica," Atti Soc. Ital. Progr. Sci. (1933) — the unilateral-contact formulation behind the active-set method.
analysis: "contact". Give the solver a list of rigid stops. Then read the contact patch that it finds. Read the docs →