verification: Coulomb friction between sliding blocks

An oracle case of a block pressed on another and slid with static and kinetic coefficients
(mu_k N carried once slipping), and the docs of the friction formulation; ROADMAP M7 and
parity updated.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01UhhZMBTviDQxZvwnyP89fo

docs/fem.md

@@ -501,16 +501,35 @@ With nonlinear geometry the pairs are searched (bucket grid of the main faces' b
 their depth limit) and integrated on the current configuration at every evaluation; a pair deeper
 than half the main element's thickness has passed through and is left out.
 
+**Friction**: Coulomb's law at the closed secondary nodes by a return map (Wriggers 2006, §5.2;
+Simo & Laursen 1992). The node's weighted tangential slip since the last converged increment is
+w_A = B_A Δx with B_A = [∫ Φ_A N_i P, −∫ Φ_A N_j P] (P = I − n n; the rate form of the mortar slip,
+McDevitt & Laursen 2000), ĝ_T = w_A / A_A; the trial traction t* = P t_n − ε_A ĝ_T sticks when
+|t*| ≤ μ p_A and otherwise returns to the cone, t = μ p_A t*/|t*|, with the static coefficient for
+a node that stuck in the last increment and the kinetic one for one that slipped (a static step
+has no slip rate for the exponential decay). The traction does the virtual work t · B_A δx, so a
+slipping interface carries exactly μ times its normal force. The tangent is symmetrised: it
+keeps (ε_A/A_A) Bᵀ B in stick and (ε_A μ p_A/(A_A |t*|)) Bᵀ (I − t̂ t̂) B in slip and leaves out
+the coupling to the pressure (μ t̂ ⊗ ∂p_A/∂x) and the change of B with the geometry, so the
+factorization stays the symmetric one, like the follower pressure's; the residual is exact, so
+only the rate of convergence suffers while nodes slip. Tractions and slip states are committed
+when an increment converges. If an increment's active set chatters (a node touching with no
+pressure, at the edge of a contact that lifts under friction), the increment is cut back and the
+nodes open by less than 10⁻³ of their element's thickness keep their penalty in the tangent for
+the rest of the step.
+
 **Initial overclosure**: nodes whose gap starts negative are reported; `adjust` takes the overlap
 out of their weighted gap (an offset: strain-free, nodes are not moved), `refuse` stops naming
 the node and the depth.
 
 **Output**: per secondary node the pressure (its normal contact force over the area it carries,
 ∫ Φ_A over the overlaps, which is smoother than the multipliers at the contact's edge), the
-multiplier, the gap ĝ and the status (open, closed); per contact the force on the secondary
+multiplier, the gap ĝ, the status (open, closed; with friction stick or slip), and with friction
+the shear traction and the slip over the last increment; per contact the force on the secondary
 surface, the normal force Σ p_A ∫Φ_A, the deepest penetration and the augmentations.
 `solve.nonlinear` summarises each contact, `solve.probe` gives a secondary node's `contact`
-entries, and the results display draws `contact_pressure`, `contact_gap` and `contact_status`.
+entries, and the results display draws `contact_pressure`, `contact_gap`, `contact_status`,
+`contact_shear` and `contact_slip`.
 
 **Verification** (`crates/solver/tests/contact.rs`, `tests/oracle/nl-contact-patch`,
 `tests/oracle/nl-hertz-cylinder`): the contact patch test on non-matching hex8 and tet4 meshes
@@ -519,7 +538,10 @@ finite sliding, where the gaps the augmentation leaves enter through their geome
 Hertz's line contact, plane strain (peak pressure within 0.3%, the pressure at two nodes within
 3% and 6%, the vertical stress under the contact within 2%, the half-width bracketed); Hertz's
 sphere on a flat, a quarter model (peak pressure within 4%, contact radius bracketed; heavy);
-Newton's quadratic rate near the solution with finite sliding; initial overclosure adjusted and
+Newton's quadratic rate near the solution with finite sliding; Coulomb friction between pressed
+blocks (stick well below the cone with the forces of glued blocks, slip on it carrying μ N to 0.1%
+with small sliding and 3% with finite sliding, where the block turns a little, and μ_k N once
+slipping; `tests/oracle/nl-friction-sliding-blocks`); initial overclosure adjusted and
 refused; bit-for-bit results on one and four threads; the refusals.
 
 ## Hyperelastic materials in explicit dynamics

docs/oracle.md

@@ -2,7 +2,7 @@
 
 Generated by `cargo xtask oracle` from [`tests/oracle/`](../tests/oracle/); the policy is [verification.md](verification.md). Every case is our own model, built with FieldCraft commands from the problem's published description, and compared with the reference value its `expected.json` cites, within the stated tolerance and for the stated reason. Accepted results are frozen in each case's `fixture.json`: a later change that moves one beyond its band (10⁻⁹ of its size) fails CI even while the reference still passes. No reference value was produced with a commercial product.
 
-**222 of 228 quantities pass** (6 known failure(s)) in 64 active cases: 55 theory, 9 published benchmark, 0 experiment; 10 pending cases wait for capabilities still being built.
+**225 of 231 quantities pass** (6 known failure(s)) in 65 active cases: 56 theory, 9 published benchmark, 0 experiment; 10 pending cases wait for capabilities still being built.
 
 ## Theory: closed-form solutions, patch tests and mesh invariants
 
@@ -105,6 +105,9 @@ Generated by `cargo xtask oracle` from [`tests/oracle/`](../tests/oracle/); the
 | [nl-elastica-cantilever](../tests/oracle/nl-elastica-cantilever/) | tip shortening δ_h/L at P L²/(E I) = 5 | 0.38763 | 0.38955 | 0.495% | 2% | pass | matches |
 | [nl-elastica-cantilever](../tests/oracle/nl-elastica-cantilever/) | tip deflection δ_v/L at P L²/(E I) = 10 | 0.81061 | 0.81421 | 0.444% | 1% | pass | matches |
 | [nl-elastica-cantilever](../tests/oracle/nl-elastica-cantilever/) | tip shortening δ_h/L at P L²/(E I) = 10 | 0.555 | 0.558 | 0.541% | 2% | pass | matches |
+| [nl-friction-sliding-blocks](../tests/oracle/nl-friction-sliding-blocks/) | the tangential over the normal force across the sliding interface | 0.2 | 0.2 | -0.000103% | 0.1% | pass | matches |
+| [nl-friction-sliding-blocks](../tests/oracle/nl-friction-sliding-blocks/) | the tangential contact force on the upper block over the sideways reaction at its top (they balance) | -1 | -1 | 0.00000224% | 1e-6 rel. | pass | matches |
+| [nl-friction-sliding-blocks](../tests/oracle/nl-friction-sliding-blocks/) | the secondary nodes that slip (all 9) | 9 | 9 | ≈ 0 | exact | pass | matches |
 | [nl-hertz-cylinder](../tests/oracle/nl-hertz-cylinder/) | the largest contact pressure p₀ (at x = 0) | 2884.6 | 2892 | 0.254% | 2% | pass | matches |
 | [nl-hertz-cylinder](../tests/oracle/nl-hertz-cylinder/) | the force across the contact (half model, 1 mm thick) | 1132.8 | 1133 | 0.0221% | 0.2% | pass | matches |
 | [nl-hertz-cylinder](../tests/oracle/nl-hertz-cylinder/) | the contact pressure at the secondary node at x = 0.29594 (node 7) | 2325.1 | 2257.8 | -2.89% | 5% | pass | matches |
@@ -283,6 +286,7 @@ Generated by `cargo xtask oracle` from [`tests/oracle/`](../tests/oracle/); the
 - **[nl-contact-patch](../tests/oracle/nl-contact-patch/)** (theory): Contact patch test: two blocks with non-matching meshes in uniaxial compression (implicit mortar contact). Source: The contact patch test: a contact discretisation must transmit a uniform pressure across a non-matching interface exactly; Taylor & Papadopoulos, On a patch test for contact problems in two dimensions, Nonlinear Computational Mechanics, Springer 1991, 690–702; Crisfield, Re-visiting the contact patch test, IJNME 48:435–449, 2000; passed by mortar methods: Puso & Laursen, CMAME 193:601–629, 2004. Uniaxial stress between frictionless symmetry planes is uniform in both blocks. Tolerances — small_lower_corner: small sliding is linear and its mortar integrals exact on these rectangular faces: a uniform stress is the discrete solution up to the convergence tolerance (10⁻¹⁰); small_lower_centre: as for the corner element; small_lower_base: as for the corner element; small_upper_top: as for the corner element; small_lateral: frictionless symmetry planes and a frictionless interface leave the blocks free to expand sideways; the non-matching interface must not hold them; small_pressure: the mortar multipliers are the nodal values of the contact pressure, uniform and equal to the stress across the interface; small_force: the pressure carries the whole load across the interface; finite_lower_corner: finite sliding is exact up to the geometric terms of the gaps the augmentation leaves (10⁻⁵ of the element thickness) and the default convergence tolerances; finite_lower_centre: as for the corner element; finite_lower_base: as for the corner element; finite_upper_top: as for the corner element; finite_lateral: frictionless symmetry planes and a frictionless interface leave the blocks free to expand sideways; the non-matching interface must not hold them; finite_pressure: the mortar multipliers are the nodal values of the contact pressure, uniform and equal to the stress across the interface; finite_force: the pressure carries the whole load across the interface; small_stress: the squeeze less the penetration the augmented Lagrangian leaves (up to 10⁻⁵ of the element thickness, a few 10⁻⁵ of the squeeze)
 - **[nl-cylinder-inflation](../tests/oracle/nl-cylinder-inflation/)** (theory): Inflation of a thick neo-Hookean cylinder by a follower pressure (plane strain). Source: Finite inflation of an incompressible neo-Hookean thick-walled tube in plane strain: r² = R² + a² − A², σ_θ − σ_r = μ (λ² − λ⁻²) with λ = r/R, integrated through the wall: p = (μ/2)[(ln λ_a² − λ_a⁻²) − (ln λ_b² − λ_b⁻²)] (R. S. Rivlin, Large elastic deformations of isotropic materials VI, Phil. Trans. R. Soc. A 242:173–195, 1949; R. W. Ogden, Non-Linear Elastic Deformations, Dover 1997, §5.3.2). Tolerances — ua: the closed form is incompressible and our material nearly so (κ/μ = 1000, a volume change of order p/κ ≈ 4·10⁻⁴); hybrid hex20, 6 × 16, follow the wall's strain gradient well inside; ub: as at the bore
 - **[nl-elastica-cantilever](../tests/oracle/nl-elastica-cantilever/)** (theory): Large-deflection cantilever under a dead tip load against the elastica (solids). Source: The elastica of a cantilever under a vertical end load: K. E. Bisshopp, D. C. Drucker, Large deflection of cantilever beams, Quarterly of Applied Mathematics 3(3):272–275, 1945 (elliptic-integral solution); values computed by integrating θ″ = −(P/EI) cos θ by shooting (our own RK4, 20000 steps), agreeing with the table of K. Mattiasson, Numerical results from large deflection beam and frame problems analysed by means of elliptic integrals, IJNME 17(1):145–153, 1981, to all six digits. Tolerances — v1: hex20, 40 × 2 through the depth, L/h = 20: shear deformation (0.6 (h/L)² ≈ 0.15 % of the deflection) and the clamp's local restraint, both absent from the elastica; u1: as δ_v; the shortening is the smaller, second-order quantity, so its relative error is about twice the deflection's; v2: hex20, 40 × 2 through the depth, L/h = 20: shear deformation (0.6 (h/L)² ≈ 0.15 % of the deflection) and the clamp's local restraint, both absent from the elastica; u2: as δ_v; the shortening is the smaller, second-order quantity, so its relative error is about twice the deflection's; v5: hex20, 40 × 2 through the depth, L/h = 20: shear deformation (0.6 (h/L)² ≈ 0.15 % of the deflection) and the clamp's local restraint, both absent from the elastica; u5: as δ_v; the shortening is the smaller, second-order quantity, so its relative error is about twice the deflection's; v10: hex20, 40 × 2 through the depth, L/h = 20: shear deformation (0.6 (h/L)² ≈ 0.15 % of the deflection) and the clamp's local restraint, both absent from the elastica; u10: as δ_v; the shortening is the smaller, second-order quantity, so its relative error is about twice the deflection's
+- **[nl-friction-sliding-blocks](../tests/oracle/nl-friction-sliding-blocks/)** (theory): Coulomb friction: a block pressed on another and slid (implicit mortar contact). Source: Coulomb's law of dry friction: a sliding interface carries μ_k times the normal force against the slip (C. A. Coulomb, Théorie des machines simples, Mém. Math. Phys. Acad. Sci. 10:161–332, 1785; P. Wriggers, Computational Contact Mechanics, 2nd ed., Springer 2006, §5.2). Tolerances — kinetic: once the interface slips its friction is μ_k N exactly (the return map puts every slipping node's traction on the cone); the force opposes the upper block's motion; balance: the friction carries the whole sideways load; slipping: a sideways demand far beyond μ N slips the whole interface
 - **[nl-hertz-cylinder](../tests/oracle/nl-hertz-cylinder/)** (theory): Hertz line contact: a cylinder pressed on a flat (implicit mortar contact, plane strain). Source: H. Hertz, Über die Berührung fester elastischer Körper, J. reine angew. Math. 92:156–171, 1882; K. L. Johnson, Contact Mechanics, Cambridge University Press 1985, §4.3 (two-dimensional contact of cylindrical bodies, eqs. 4.45–4.49): b = √(4PR/(πE*)), p₀ = 2P/(πb), with 1/E* = 2(1 − ν²)/E for two bodies of one material and R the cylinder's radius on a flat. Tolerances — peak_pressure: p₀ = 2P/(πb) for the load per length P = 2265.57 N/mm, b = 0.5 mm; the contact width spans seven elements of the graded mesh; normal_force: the pressure carries the whole load; the soft springs that hold the indenter before contact take a millionth; pressure_x0296: p₀ √(1 − x²/b²) at the node; the nodal pressure (the node's normal force over the area it carries) averages the distribution over the node's elements; pressure_x0367: as at node 7, nearer the edge, where the elements round off the square-root fall of the pressure to zero at b; pressure_x0633: outside the contact the surfaces are apart: no pressure (an element beyond b, since the edge element carries part of the last closed node's pressure); syy_depth: σ_z = −p₀ / √(1 + z²/b²) on the axis of a half-plane (Johnson 1985, eq. 4.49); the foundation is ten half-widths wide and deep, not a half-plane, and the stress is the nodal average of the elements around the node
 - **[nl-hybrid-incompressible](../tests/oracle/nl-hybrid-incompressible/)** (theory): Nearly incompressible hybrid block under prescribed stretch. Source: Uniaxial stress of the decoupled compressible neo-Hookean energy W = μ/2 (Ī₁ − 3) + κ/2 (J − 1)² (P. J. Flory, Trans. Faraday Soc. 57:829–838, 1961): σ = μ J^(−5/3) (λ² − I₁/3) + κ (J − 1) with the lateral stretch found so σ_yy = 0 (bisection; E = 3, ν = 0.4999: μ = 1.0000667, κ = 5000). Tolerances — sigma: a homogeneous deformation every hybrid hex8 represents exactly; the 27-element mesh makes the global mixed tangent indefinite in the first iterations (LDLᵀ): the tolerance admits the Newton tolerance and rounding; lateral: a homogeneous deformation every hybrid hex8 represents exactly; the 27-element mesh makes the global mixed tangent indefinite in the first iterations (LDLᵀ): the tolerance admits the Newton tolerance and rounding
 - **[nl-j2-uniaxial](../tests/oracle/nl-j2-uniaxial/)** (theory): Uniaxial tension of a J2 bilinear-hardening cube to large strain (implicit, finite strain). Source: J2 plasticity with linear isotropic hardening in uniaxial stress (Simo and Hughes, Computational Inelasticity, Springer 1998, §2.3) in the logarithmic strain, which Simo's exponential return map (CMAME 99:61–112, 1992) integrates exactly on a proportional path: Kirchhoff stress τ = σ_y + E_t (ε − σ_y/E), Cauchy stress σ = τ/J with ln J = τ/(3K) (Hencky's volume law). Tolerances — sigma: a homogeneous deformation the element represents exactly, reached on a proportional path: the tolerance only admits the Newton tolerance and rounding; plastic_strain: ε_p = (τ − σ_y)/H, as for the stress

docs/parity.md

@@ -55,7 +55,7 @@ products are run, no documentation text is copied). Status: **done**, **started*
 | Model | Steps: linear static, frequency; output presets | M3 | | done |
 | Model | Steps: nonlinear static, explicit dynamics, heat transfer; amplitudes | M7–M8, M15 | | started (nonlinear static with or without nonlinear geometry, increments, tolerances and a load path (cyclic loading), written to Abaqus `*STEP, NLGEOM` and LS-DYNA implicit (the load path as LS-DYNA's curve; refused by Abaqus for now), refused by Nastran; explicit dynamics; heat transfer; amplitudes per load planned) |
 | Interaction | Tie between non-matching meshes | M7 | | planned |
-| Interaction | Contact (general and pairs), friction, automatic pair detection | M7, M10 | | started (explicit: pairs, general and self contact, segment to segment with a penalty that keeps the time increment, Coulomb friction with exponential decay, friction on rigid walls, `interaction.contact.find`, Abaqus and LS-DYNA contact keywords; implicit: frictionless mortar contact pairs in nonlinear static steps by an augmented Lagrangian, friction and self contact next) |
+| Interaction | Contact (general and pairs), friction, automatic pair detection | M7, M10 | | started (explicit: pairs, general and self contact, segment to segment with a penalty that keeps the time increment, Coulomb friction with exponential decay, friction on rigid walls, `interaction.contact.find`, Abaqus and LS-DYNA contact keywords; implicit: mortar contact pairs in nonlinear static steps by an augmented Lagrangian with Coulomb friction, self contact next) |
 | Interaction | Springs and dashpots (ground or two-point; DOF, vector or line), point masses with rotary inertia | M4 | | done (implicit; explicit translational; all three decks with what each cannot take reported) |
 | Interaction | Bolt pretension, connectors | M7 | | planned |
 | Decks | Abaqus input writer (`deck.write`, flat deck; reports what Abaqus cannot express) | M3 | | done |
@@ -63,7 +63,7 @@ products are run, no documentation text is copied). Status: **done**, **started*
 | Decks | LS-DYNA keyword writer (implicit; explicit in M8) | M3 | | done (implicit static and eigenvalue; beam sections planned) |
 | Decks | Readers keeping unknown keywords | M9 | | planned |
 | Solve | Own linear static and modal solver, verified against analytic solutions | M4 | yes | started (solid elements: tet4, tet10, wedge6, hex8, hex20; 2D plane stress, plane strain and axisymmetric tri3, tri6, quad4 (full, incompatible modes), quad8 (full, reduced); quad4/tri3 MITC shells and Timoshenko beams with profiles, mixed with solids; reference points, kinematic and distributing couplings, moments; centrifugal loads) |
-| Solve | Own nonlinear implicit solver: Newton–Raphson with line search and automatic increments, total Lagrangian solids and 2D elements, mixed u–p at finite strain, follower pressure | M15 | yes | started (solids, plane strain, axisymmetric; St Venant–Kirchhoff, neo-Hookean and J2 plasticity (Simo's exponential return map at finite strain, radial return at small strain, consistent tangents, the return shared with explicit dynamics); verified by the elastica, tube inflation, patch and objectivity tests, uniaxial tension to large strain, Hill's thick cylinder, a cantilever's collapse and the Bauschinger effect; frictionless mortar contact (patch test, Hertz); Mooney–Rivlin and Yeoh; thermal strain at finite and small strain; a symmetric indefinite LDLᵀ factorization with inertia, Newton through indefinite tangents (hybrid elements at ν = 0.4999 under prescribed displacement); linear buckling of solids (Euler columns, Bryan's plate) and imperfections from its modes; arc-length with critical points from the inertia (shallow-arch snap-through); shells and implicit friction next) |
+| Solve | Own nonlinear implicit solver: Newton–Raphson with line search and automatic increments, total Lagrangian solids and 2D elements, mixed u–p at finite strain, follower pressure | M15 | yes | started (solids, plane strain, axisymmetric; St Venant–Kirchhoff, neo-Hookean and J2 plasticity (Simo's exponential return map at finite strain, radial return at small strain, consistent tangents, the return shared with explicit dynamics); verified by the elastica, tube inflation, patch and objectivity tests, uniaxial tension to large strain, Hill's thick cylinder, a cantilever's collapse and the Bauschinger effect; mortar contact with friction (patch test, Hertz, sliding blocks); Mooney–Rivlin and Yeoh; thermal strain at finite and small strain; a symmetric indefinite LDLᵀ factorization with inertia, Newton through indefinite tangents (hybrid elements at ν = 0.4999 under prescribed displacement); linear buckling of solids (Euler columns, Bryan's plate) and imperfections from its modes; arc-length with critical points from the inertia (shallow-arch snap-through); shells next) |
 | Solve | Own explicit dynamics solver (central differences, one-point hex8 with hourglass control, nodal-pressure tet4, bulk viscosity, energy balance), verified against analytic wave and impact solutions | M10 | yes | started (X1 core and X2: large deformation, objective rates, J2 and Johnson–Cook plasticity, damage, erosion; X3: deformable contact and friction, verified by bar collision, sliding blocks, the contact patch test and Hertz; hyperelastic materials on the shared material point; opt-in mass scaling) |
 | Results | Deformed shape, contours, mode shapes, probes, reaction sums | M4 | | started (results in the session, `solve.probe`, reaction sums; display planned) |
 | Results | XY and path plots, cuts, display groups, free-body forces, animation | M9 | | planned |

docs/solvers/literature.md

@@ -1322,6 +1322,13 @@ written from (cited again in the module docs); docs/fem.md describes the formula
   Laursen, *An augmented Lagrangian treatment of contact problems involving friction*, Comput.
   Struct. 42(1):97–116, 1992; P. Wriggers, *Computational Contact Mechanics*, 2nd ed., Springer
   2006, ch. 6 and 10; T. A. Laursen, *Computational Contact and Impact Mechanics*, Springer 2003.
+- Friction: the regularised Coulomb return map at the nodes, J. C. Simo, T. A. Laursen 1992 and
+  P. Wriggers 2006, §5.2 (as in Part Q), on the mortar slip in its rate form (the weighted
+  relative motion projected on the tangent plane), T. W. McDevitt, T. A. Laursen 2000 and M. A.
+  Puso, T. A. Laursen, *A mortar segment-to-segment frictional contact method for large
+  deformations*, CMAME 193:4891–4913, 2004; C. A. Coulomb's law, *Théorie des machines simples*,
+  1785. The tangent is symmetrised like the follower pressure's (the non-symmetric coupling to
+  the pressure is left out; the residual is exact).
 - The active set inside Newton (the max() of the regularised pressure, a semismooth Newton
   method): P. Alart, A. Curnier, *A mixed formulation for frictional contact problems prone to
   Newton like solution methods*, CMAME 92:353–375, 1991; Popp et al. 2009.