From 6ca3233aaaa9eb4cacd309fd7191cd689afc5953 Mon Sep 17 00:00:00 2001 From: jyviko <2501648+jyviko@users.noreply.github.com> Date: Wed, 14 Oct 2026 15:14:04 -0400 Subject: [PATCH] solver: Coulomb friction for implicit mortar contact A return map at the closed secondary nodes on the mortar slip in its rate form (w_A = B_A dx, B_A = [int Phi_A N_i P, -int Phi_A N_j P]): stick within mu p_A, slip back to the cone, static coefficient for a node that stuck in the last increment and kinetic for one that slipped; the friction does the virtual work t . B_A dx, so a slipping interface carries mu N exactly. The tangent is symmetrised (stick and slip stiffness kept, the coupling to the pressure and the change of B left out; the residual is exact). Tractions and slip states are committed with the converged increment. After a chattering cutback, nodes open by less than 1e-3 of their element's thickness keep their penalty in the tangent for the rest of the step. Results carry the shear traction, the slip and stick/slip status. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01UhhZMBTviDQxZvwnyP89fo --- crates/model/src/contact.rs | 2 + crates/solver/src/contact.rs | 370 +++++++++++++++++++++++++++++++-- crates/solver/src/nonlinear.rs | 66 ++++-- crates/solver/tests/contact.rs | 60 ++++++ 4 files changed, 463 insertions(+), 35 deletions(-) diff --git a/crates/model/src/contact.rs b/crates/model/src/contact.rs index 0f992c1..ce8a051 100644 --- a/crates/model/src/contact.rs +++ b/crates/model/src/contact.rs @@ -13,6 +13,8 @@ //! **Friction** is Coulomb's law with a coefficient that falls from its static to its kinetic //! value with the slip rate γ̇: μ = μ_k + (μ_s − μ_k) e^{−d γ̇} (the exponential decay of //! E. Rabinowicz's observations, as LS-DYNA's and Abaqus's exponential decay friction write it). +//! A static step has no slip rate: its nodes take the static coefficient while they stick and the +//! kinetic one once they slip. //! The same friction applies to the planar rigid walls of an explicit step //! ([`crate::RigidWall`]). //! diff --git a/crates/solver/src/contact.rs b/crates/solver/src/contact.rs index d6678f0..41e7bf1 100644 --- a/crates/solver/src/contact.rs +++ b/crates/solver/src/contact.rs @@ -76,7 +76,7 @@ //! per contact the force on the secondary surface and its normal part. //! //! **Scope.** Pairs of surfaces of linear solid faces (tri3 and quad4 faces: tet4, wedge6 and -//! hex8 elements) in 3D models, frictionless. General (single-surface) contact, quadratic faces, +//! hex8 elements) in 3D models, with Coulomb friction ([`Interface::friction`]). General (single-surface) contact, quadratic faces, //! 2D models and arc-length steps are refused by name. use std::collections::BTreeMap; @@ -85,7 +85,7 @@ use fieldcraft_fem::contact::{Aabb, Grid, bounds}; use fieldcraft_fem::dual::{Dual, Hyper, Real, VARS}; use fieldcraft_fem::frame::{Vec3, dot, norm, unit}; use fieldcraft_fem::mortar::{area_vector, integrate}; -use fieldcraft_model::{Contact, ContactSurfaces, ElementKind, InitialPenetration, Model, Resolved}; +use fieldcraft_model::{Contact, ContactSurfaces, ElementKind, Friction, InitialPenetration, Model, Resolved}; use fieldcraft_parallel as parallel; use serde::Serialize; @@ -105,10 +105,19 @@ pub const PENALTY_FACTOR: f64 = 10.0; /// the element under it (and no open node keeps a multiplier). pub const GAP_TOLERANCE: f64 = 1e-5; +/// Once an increment's active set has chattered, a node open by less than this fraction of its +/// element's thickness keeps its penalty in the tangent for the rest of the step (its +/// pressure, and so the residual, is still zero): a node that touches with no pressure, the edge +/// of a contact that lifts, otherwise flips between open and closed from one Newton iteration to +/// the next, each state's correction pushing it into the other, and a smaller increment does not +/// help under proportional loading. +pub const ENGAGED_BAND: f64 = 1e-3; + /// Most augmentations of one increment. pub const MAX_AUGMENTATIONS: usize = 25; -/// An increment is cut back when the active set alternates between two states this many times. +/// An increment is cut back when the active set alternates between two states this many times +/// (and the step takes [`ENGAGED_BAND`] from then on). pub const CHATTER_LIMIT: usize = 3; /// Most faces over the contacts of a model. @@ -149,20 +158,48 @@ pub struct Interface { /// fixed in the reference configuration; `None`: finite sliding, searched and integrated /// on the current configuration at every evaluation. fixed: Option>>, + /// Coulomb friction (frictionless: none). + friction: Option, } -/// Whether a secondary node presses on the main surface. +/// A contact's state carried from one converged increment to the next: per secondary node the +/// augmented Lagrange multiplier, and with friction the tangential traction and whether the node +/// was slipping. +#[derive(Debug, Clone, Default, PartialEq)] +pub struct Committed { + pub multipliers: Vec, + pub traction: Vec, + pub slipping: Vec, +} + +impl Committed { + /// Take a converged state's friction (the multipliers are the augmentation's). + pub fn commit(&mut self, state: &State) { + self.traction.clone_from(&state.traction); + self.slipping.clone_from(&state.slipping); + } +} + +/// Whether a secondary node presses on the main surface (and with friction, sticks or slips). #[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize)] #[serde(rename_all = "snake_case")] pub enum Status { Open, Closed, + Stick, + Slip, +} + +impl Status { + pub fn is_closed(self) -> bool { + self != Status::Open + } } impl ContactResult { /// The closed nodes. pub fn closed(&self) -> usize { - self.status.iter().filter(|s| **s == Status::Closed).count() + self.status.iter().filter(|s| s.is_closed()).count() } } @@ -178,6 +215,12 @@ pub struct State { /// (its faces' area vectors summed). pub nodal_force: Vec, pub normal: Vec, + /// Friction: the tangential traction on each secondary node (on it, against its slip), its + /// slip over the increment (relative to the main surface, normalised like the gap) and + /// whether it slips; empty without friction. + pub traction: Vec, + pub slip: Vec, + pub slipping: Vec, } impl State { @@ -206,6 +249,12 @@ pub struct ContactResult { pub pressure: Vec, /// The mortar multipliers (the coefficients of the pressure field, Φ_A-interpolated). pub multiplier: Vec, + /// Friction: the shear traction (the tangential contact force over the area, on the node) + /// and the slip over the last increment, per node; empty without friction. + #[serde(skip_serializing_if = "Vec::is_empty")] + pub shear: Vec, + #[serde(skip_serializing_if = "Vec::is_empty")] + pub slip: Vec, /// The normalised gap ĝ (negative: penetration), `None` where no main face is over the node. pub gap: Vec>, pub status: Vec, @@ -301,6 +350,23 @@ struct Fixed { } impl Fixed { + /// B_a, the slip operator ([`slip_operator`]) of the reference configuration. + fn slip_jacobian(&self, a: usize) -> [[f64; VARS]; 3] { + let mut j = [[0.0; VARS]; 3]; + let (Some(d), Some(m)) = (self.d.get(a), self.m.get(a)) else { return j }; + let n = self.normal; + for k in 0..3 { + for l in 0..3 { + let proj = if k == l { 1.0 } else { 0.0 } - n[k] * n[l]; + for i in 0..4 { + j[k][3 * i + l] = d[i] * proj; + j[k][12 + 3 * i + l] = -m[i] * proj; + } + } + } + j + } + /// ∂g̃_a/∂x over the pair's variables (secondary corners from 0, main corners from 12). fn gradient(&self, a: usize) -> [f64; VARS] { let mut g = [0.0; VARS]; @@ -335,6 +401,35 @@ fn weighted_gaps(xs: &[[T; 3]], xm: &[[T; 3]]) -> Option<[T; 4]> { Some(g) } +/// The slip operator of one face pair (corners `xs`, `xm`): per secondary corner a the matrix +/// B_a with w_a = B_a Δx, the weighted tangential slip ∫ Φ_a (I − n n)(Δx_s − Δx_m) dA of the +/// pair's corners' motions Δx (rows: its components; columns: the pair's variables, secondary +/// corners from 0, main corners from 12), with the integrals on the current configuration: the +/// rate form of the mortar slip, B_a = [∫ Φ_a N_i P, −∫ Φ_a N_j P] (McDevitt & Laursen 2000; Puso & +/// Laursen 2004), P = I − n n. +fn slip_operator(xs: &[Vec3], xm: &[Vec3]) -> Option<[[[f64; VARS]; 3]; 4]> { + let mut b = [[[0.0; VARS]; 3]; 4]; + integrate(xs, xm, |p| { + let n = p.normal; + for a in 0..4 { + let wa = p.weight * p.secondary[a]; + if wa == 0.0 { + continue; + } + for k in 0..3 { + for l in 0..3 { + let proj = wa * (if k == l { 1.0 } else { 0.0 } - n[k] * n[l]); + for c in 0..4 { + b[a][k][3 * c + l] += proj * p.secondary[c]; + b[a][k][12 + 3 * c + l] -= proj * p.main[c]; + } + } + } + } + })?; + Some(b) +} + /// The variables of a face pair: secondary corners at 0, main corners from 12. fn variables(xs: &[Vec3], xm: &[Vec3], var: impl Fn(f64, usize) -> T) -> (Vec<[T; 3]>, Vec<[T; 3]>) { let s = xs.iter().enumerate().map(|(i, x)| [0, 1, 2].map(|k| var(x[k], 3 * i + k))).collect(); @@ -368,9 +463,6 @@ impl Interface { "general (single-surface) contact `{name}` in implicit steps; define contact pairs (interaction.contact.pair), or solve it in explicit dynamics" ))); }; - if !contact.friction.is_frictionless() { - return Err(Error::NotSupported(format!("friction (contact `{name}`) in implicit steps"))); - } let secondary = faces(model, resolved, mesh, name, secondary)?; let main = faces(model, resolved, mesh, name, main)?; let mut nodes: Vec = secondary.iter().flat_map(|f| f.nodes.iter().copied()).collect(); @@ -406,7 +498,18 @@ impl Interface { return Err(Error::Invalid(format!("contact `{name}`: the secondary faces at node {label} have no area"))); } let offset = vec![0.0; nodes.len()]; - Ok(Interface { name: name.to_string(), secondary, main, nodes, area, penalty, thickness, offset, fixed: None }) + Ok(Interface { + name: name.to_string(), + secondary, + main, + nodes, + area, + penalty, + thickness, + offset, + fixed: None, + friction: (!contact.friction.is_frictionless()).then_some(contact.friction), + }) } /// Adjust (offset) or refuse the nodes that start behind the main surface. @@ -577,13 +680,47 @@ impl Interface { (gap, support, index) } + /// The state a step starts from: no multipliers, no friction traction. + pub fn committed(&self) -> Committed { + let n = self.nodes.len(); + let friction = self.friction.is_some(); + Committed { + multipliers: vec![0.0; n], + traction: if friction { vec![[0.0; 3]; n] } else { Vec::new() }, + slipping: vec![false; if friction { n } else { 0 }], + } + } + /// The contact's forces at displacements `u` (added to `internal`), its tangent blocks when - /// `tangent`, and its state, with the multipliers `multipliers` (per node). - pub fn evaluate(&self, mesh: &Mesh, u: &[f64], multipliers: &[f64], tangent: bool, internal: &mut [f64], blocks: &mut Blocks) -> Result { + /// `tangent`, and its state, from the last converged increment's state `committed` and + /// displacements `reference`; nodes open by less than `band` of their element's thickness + /// keep their penalty in the tangent ([`ENGAGED_BAND`]; 0: only closed nodes). + #[allow(clippy::too_many_arguments)] + pub fn evaluate( + &self, + mesh: &Mesh, + u: &[f64], + reference: &[f64], + committed: &Committed, + band: f64, + tangent: bool, + internal: &mut [f64], + blocks: &mut Blocks, + ) -> Result { + let multipliers = committed.multipliers.as_slice(); let pairs = self.pairs(mesh, u)?; let (weighted, support, corner_index) = self.sums(&pairs); let n = self.nodes.len(); - let mut state = State { gap: vec![None; n], pressure: vec![0.0; n], support, nodal_force: vec![[0.0; 3]; n], normal: vec![[0.0; 3]; n] }; + let mut state = State { + gap: vec![None; n], + pressure: vec![0.0; n], + support, + nodal_force: vec![[0.0; 3]; n], + normal: vec![[0.0; 3]; n], + traction: Vec::new(), + slip: Vec::new(), + slipping: Vec::new(), + }; for (f, corners) in self.secondary.iter().zip(&corner_index) { let Some(a) = area_vector(&Interface::corners(mesh, u, f)) else { continue }; for i in corners { @@ -595,19 +732,23 @@ impl Interface { for v in &mut state.normal { *v = unit(*v).unwrap_or(*v); } + // Nodes closed, or open by less than `band` of their element's thickness: their penalty + // stays in the tangent. + let mut engaged = vec![false; n]; for i in 0..self.nodes.len() { let (Some(Some(g)), Some(a), Some(e)) = (weighted.get(i), self.area.get(i), self.penalty.get(i)) else { continue }; let gap = g / a; let lambda = multipliers.get(i).copied().unwrap_or(0.0); - if let (Some(slot), Some(p)) = (state.gap.get_mut(i), state.pressure.get_mut(i)) { + if let (Some(slot), Some(p), Some(h)) = (state.gap.get_mut(i), state.pressure.get_mut(i), self.thickness.get(i)) { *slot = Some(gap); *p = (lambda - e * gap).max(0.0); + engaged[i] = lambda - e * gap > -e * band * h; } } // The pairs with a closed secondary corner, linearised. let mut active: Vec<(usize, usize)> = Vec::new(); for (s, found) in pairs.iter().enumerate() { - let closed = corner_index.get(s).is_some_and(|c| c.iter().any(|i| state.pressure.get(*i).is_some_and(|p| *p > 0.0))); + let closed = corner_index.get(s).is_some_and(|c| c.iter().any(|i| engaged.get(*i).copied().unwrap_or(false))); if closed { active.extend((0..found.len()).map(|k| (s, k))); } @@ -663,7 +804,7 @@ impl Interface { if let Some(own) = own { (0..3).for_each(|k| own[k] += on[k]); } - if p > 0.0 && tangent { + if tangent && engaged.get(*i).copied().unwrap_or(false) { let entry = node_gradients.entry(*i).or_default().entry(*node).or_insert([0.0; 3]); for k in 0..3 { entry[k] += gv[k]; @@ -710,22 +851,213 @@ impl Interface { } blocks.push((nodes, data)); } + if let Some(law) = self.friction { + self.friction(mesh, u, reference, committed, (&pairs, &corner_index, &active), law, tangent, &mut state, internal, blocks)?; + } Ok(state) } + /// Coulomb friction at the closed secondary nodes by a return map (Wriggers 2006, §5.2; + /// Simo & Laursen 1992): the weighted tangential slip of each node since the last converged + /// increment, w_A = Σ ∫ Φ_A (I − n n)(Δx_s − Δx_m) dA = B_A Δx ([`slip_operator`]) and + /// ĝ_T = w_A / A_A; the trial traction + /// t* = P t_n − ε_A ĝ_T (P projects the committed traction on the current tangent plane); + /// stick, t = t*, when |t*| ≤ μ p_A, else slip back to the cone, t = μ p_A t*/|t*|, μ the + /// static coefficient for a node that stuck in the last increment and the kinetic one for + /// one that slipped. The traction does the virtual work t · B_A δx, so the forces −B_Aᵀ t enter + /// the internal forces (B_A at the current configuration: the friction force is exactly μ + /// times the normal one where nodes slip). The tangent keeps the symmetric stick/slip + /// stiffness (ε_A/A_A) Bᵀ B, or (ε_A μ p_A/(A_A |t*|)) Bᵀ (I − t̂ t̂) B for slip, and leaves out + /// the coupling to the pressure (μ t̂ ⊗ ∂p_A/∂x, not symmetric) and the change of B with the + /// geometry: the factorization stays symmetric and the residual exact, so only the rate of + /// convergence can suffer while nodes slip. + #[allow(clippy::too_many_arguments)] + fn friction( + &self, + mesh: &Mesh, + u: &[f64], + reference: &[f64], + committed: &Committed, + (pairs, corner_index, active): (&[Vec], &[Vec], &[(usize, usize)]), + law: Friction, + tangent: bool, + state: &mut State, + internal: &mut [f64], + blocks: &mut Blocks, + ) -> Result<()> { + let n = self.nodes.len(); + state.traction = vec![[0.0; 3]; n]; + state.slip = vec![[0.0; 3]; n]; + state.slipping = vec![false; n]; + // Per active pair: its corners' weighted slips and their Jacobians. + let terms = parallel::map_grained(active.len(), GRAIN, |t| { + let &(s, k) = active.get(t)?; + let (f, corners) = (self.secondary.get(s)?, corner_index.get(s)?); + let m = self.main.get(pairs.get(s)?.get(k)?.main)?; + let (xs, xm) = (Interface::corners(mesh, u, f), Interface::corners(mesh, u, m)); + let (rs, rm) = (Interface::corners(mesh, reference, f), Interface::corners(mesh, reference, m)); + let operator: Vec<[[f64; VARS]; 3]> = match self.fixed.as_ref().and_then(|f| f.get(s)).and_then(|l| l.get(k)) { + Some(fixed) => (0..corners.len()).map(|a| fixed.slip_jacobian(a)).collect(), + None => slip_operator(&xs, &xm)?.into_iter().take(corners.len()).collect(), + }; + let mut delta = [0.0; VARS]; + for (i, (x, r)) in xs.iter().zip(&rs).enumerate() { + (0..3).for_each(|c| delta[3 * i + c] = x[c] - r[c]); + } + for (j, (x, r)) in xm.iter().zip(&rm).enumerate() { + (0..3).for_each(|c| delta[12 + 3 * j + c] = x[c] - r[c]); + } + let slips: Vec<(Vec3, [[f64; VARS]; 3])> = + operator.into_iter().map(|b| (b.map(|row| row.iter().zip(&delta).map(|(x, y)| x * y).sum()), b)).collect(); + let nodes: Vec = f.nodes.iter().chain(&m.nodes).copied().collect(); + let dofs: Vec<[usize; 3]> = f.dofs.iter().chain(&m.dofs).copied().collect(); + Some((nodes, dofs, corners.clone(), slips)) + }); + let terms: Vec<_> = terms.into_iter().flatten().collect(); + let var = |count: usize, c: usize, l: usize| if c < count { 3 * c + l } else { 12 + 3 * (c - count) + l }; + // Node sums: the slip, and with the tangent J_A over the nodes the node reaches. + let mut total = vec![[0.0; 3]; n]; + let mut jacobians: BTreeMap> = BTreeMap::new(); + for (nodes, _, corners, slips) in &terms { + for (a, i) in corners.iter().enumerate() { + let (Some((w, jac)), Some(sum)) = (slips.get(a), total.get_mut(*i)) else { continue }; + (0..3).for_each(|k| sum[k] += w[k]); + if tangent && state.pressure.get(*i).is_some_and(|p| *p > 0.0) { + for (c, node) in nodes.iter().enumerate() { + let block = jacobians.entry(*i).or_default().entry(*node).or_insert([[0.0; 3]; 3]); + for k in 0..3 { + for l in 0..3 { + block[k][l] += jac[k].get(var(corners.len(), c, l)).copied().unwrap_or(0.0); + } + } + } + } + } + } + // The return map, and per node the tangent's middle matrix. + let mut middle: Vec> = vec![None; n]; + for i in 0..n { + let (Some(&p), Some(&area), Some(&eps), Some(&normal)) = + (state.pressure.get(i), self.area.get(i), self.penalty.get(i), state.normal.get(i)) + else { + continue; + }; + if !(p > 0.0) { + continue; + } + let old = committed.traction.get(i).copied().unwrap_or_default(); + let along = dot(old, normal); + let slip = total.get(i).copied().unwrap_or_default().map(|w| w / area); + let trial: Vec3 = std::array::from_fn(|k| old[k] - along * normal[k] - eps * slip[k]); + let mu = if committed.slipping.get(i).copied().unwrap_or(false) { law.kinetic_coefficient() } else { law.static_coefficient }; + let (limit, size) = (mu * p, norm(trial)); + let (traction, slipping, m) = if size <= limit { + (trial, false, std::array::from_fn(|k| std::array::from_fn(|l| if k == l { eps / area } else { 0.0 }))) + } else { + let dir = trial.map(|c| c / size); + let c = eps * limit / (size * area); + ( + dir.map(|d| limit * d), + true, + std::array::from_fn(|k| std::array::from_fn(|l| c * (if k == l { 1.0 } else { 0.0 } - dir[k] * dir[l]))), + ) + }; + if let (Some(t), Some(g), Some(s), Some(mm)) = + (state.traction.get_mut(i), state.slip.get_mut(i), state.slipping.get_mut(i), middle.get_mut(i)) + { + (*t, *g, *s, *mm) = (traction, slip, slipping, Some(m)); + } + } + // Forces −t · ∂w/∂x. + for (_, dofs, corners, slips) in &terms { + for (a, i) in corners.iter().enumerate() { + let (Some(t), Some((_, jac))) = (state.traction.get(*i).copied(), slips.get(a)) else { continue }; + if t == [0.0; 3] { + continue; + } + for (c, d) in dofs.iter().enumerate() { + let on: Vec3 = std::array::from_fn(|l| (0..3).map(|k| t[k] * jac[k].get(var(corners.len(), c, l)).copied().unwrap_or(0.0)).sum()); + for l in 0..3 { + if let Some(x) = internal.get_mut(d[l]) { + *x -= on[l]; + } + } + if let Some(own) = corners.get(c).and_then(|j| state.nodal_force.get_mut(*j)) { + (0..3).for_each(|l| own[l] += on[l]); + } + } + } + } + // Jᵀ M J per closed node. + for (i, reach) in jacobians { + let Some(Some(m)) = middle.get(i) else { continue }; + let nodes: Vec = reach.keys().copied().collect(); + let blocks_of: Vec<&[[f64; 3]; 3]> = reach.values().collect(); + let width = 3 * nodes.len(); + let mut data = vec![0.0; width * width]; + for (a, ja) in blocks_of.iter().enumerate() { + for (b, jb) in blocks_of.iter().enumerate() { + for l in 0..3 { + for r in 0..3 { + let mut v = 0.0; + for k in 0..3 { + for q in 0..3 { + v += ja[k][l] * m[k][q] * jb[q][r]; + } + } + if let Some(x) = data.get_mut((3 * a + l) * width + 3 * b + r) { + *x = v; + } + } + } + } + } + blocks.push((nodes, data)); + } + Ok(()) + } + /// The result of a converged state. pub fn result(&self, mesh: &Mesh, state: &State, augmentations: usize) -> ContactResult { - let status: Vec = state.pressure.iter().map(|p| if *p > 0.0 { Status::Closed } else { Status::Open }).collect(); - let penetration = state.gap.iter().zip(&status).filter(|(_, s)| **s == Status::Closed).filter_map(|(g, _)| *g).fold(0.0f64, |m, g| m.max(-g)); - // The force on the secondary nodes pushes them back along −n. + let status: Vec = state + .pressure + .iter() + .enumerate() + .map(|(i, p)| match (*p > 0.0, self.friction.is_some(), state.slipping.get(i).copied().unwrap_or(false)) { + (false, ..) => Status::Open, + (true, false, _) => Status::Closed, + (true, true, false) => Status::Stick, + (true, true, true) => Status::Slip, + }) + .collect(); + let penetration = state.gap.iter().zip(&status).filter(|(_, s)| s.is_closed()).filter_map(|(g, _)| *g).fold(0.0f64, |m, g| m.max(-g)); + // The force on the secondary nodes pushes them back along −n; its tangential part is the + // friction's. let normal: Vec = state.nodal_force.iter().zip(&state.normal).map(|(f, n)| -dot(*f, *n)).collect(); - let pressure = normal.iter().zip(&state.support).map(|(f, a)| if *a > 0.0 { f / a } else { 0.0 }).collect(); + let per_area = |v: f64, a: f64| if a > 0.0 { v / a } else { 0.0 }; + let pressure = normal.iter().zip(&state.support).map(|(f, a)| per_area(*f, *a)).collect(); + let shear = if self.friction.is_some() { + state + .nodal_force + .iter() + .zip(&state.normal) + .zip(&state.support) + .map(|((f, n), a)| { + let along = dot(*f, *n); + std::array::from_fn(|k| per_area(f[k] - along * n[k], *a)) + }) + .collect() + } else { + Vec::new() + }; let force = state.nodal_force.iter().fold([0.0; 3], |s, f| [s[0] + f[0], s[1] + f[1], s[2] + f[2]]); ContactResult { name: self.name.clone(), nodes: self.nodes.iter().map(|n| mesh.labels.get(*n).copied().unwrap_or(0)).collect(), pressure, multiplier: state.pressure.clone(), + shear, + slip: state.slip.clone(), gap: state.gap.clone(), status, force, diff --git a/crates/solver/src/nonlinear.rs b/crates/solver/src/nonlinear.rs index fedb1f3..3ea8dcb 100644 --- a/crates/solver/src/nonlinear.rs +++ b/crates/solver/src/nonlinear.rs @@ -414,9 +414,19 @@ impl Problem<'_> { /// The forces at displacements `u` and load factor λ, and the tangent (symmetric part of the /// load stiffness included) when `tangent`. - /// `multipliers` are the contacts' augmented Lagrange multipliers per secondary node (none: - /// zero). - fn evaluate(&self, u: &[f64], lambda: f64, tangent: bool, committed: &[Vec], multipliers: &[Vec]) -> Result { + /// `contact` is each contact's state at the last converged increment (augmented Lagrange + /// multipliers, friction tractions; none: zero) and `reference` the displacements there (the + /// start of the friction's slip increment; none: zero). + fn evaluate( + &self, + u: &[f64], + lambda: f64, + tangent: bool, + committed: &[Vec], + contact: &[contact::Committed], + reference: &[f64], + band: f64, + ) -> Result { let layout = &self.mesh.layout; let n = u.len(); let mut internal = vec![0.0; n]; @@ -424,7 +434,8 @@ impl Problem<'_> { let mut blocks = contact::Blocks::new(); let mut states = Vec::with_capacity(self.contacts.len()); for (i, c) in self.contacts.iter().enumerate() { - states.push(c.evaluate(self.mesh, u, multipliers.get(i).map(Vec::as_slice).unwrap_or_default(), tangent, &mut internal, &mut blocks)?); + let empty = contact::Committed::default(); + states.push(c.evaluate(self.mesh, u, reference, contact.get(i).unwrap_or(&empty), band, tangent, &mut internal, &mut blocks)?); } let mut k = tangent.then(|| { if blocks.is_empty() { @@ -569,7 +580,12 @@ struct Run<'a> { undeformed: Option>, /// The contacts' multipliers per secondary node (the last converged increment's, or the /// current augmentation's), and the augmentations over the step. - multipliers: Vec>, + contact: Vec, + /// The displacements of the last converged increment (all degrees of freedom). + reference: Vec, + /// The band of the engaged contact nodes ([`contact::ENGAGED_BAND`] once an increment's + /// active set has chattered, for the rest of the step; else 0). + band: f64, augmentations: usize, /// The augmentation being solved in the current increment. augmentation: usize, @@ -611,7 +627,15 @@ impl Run<'_> { } fn state(&self, q: &[f64], lambda: f64, tangent: bool) -> std::result::Result { - Ok(self.problem.evaluate(&self.problem.red.expand_scaled(q, lambda), lambda, tangent, &self.committed, &self.multipliers)?) + Ok(self.problem.evaluate( + &self.problem.red.expand_scaled(q, lambda), + lambda, + tangent, + &self.committed, + &self.contact, + &self.reference, + self.band, + )?) } /// How the tangent of the undeformed, unloaded model factors (Cholesky): a first tangent @@ -621,7 +645,7 @@ impl Run<'_> { if let Some(u) = &self.undeformed { return Ok(u.clone()); } - let mut e = self.problem.evaluate(&vec![0.0; self.problem.mesh.dofs()], 0.0, true, &self.committed, &[])?; + let mut e = self.problem.evaluate(&vec![0.0; self.problem.mesh.dofs()], 0.0, true, &self.committed, &[], &[], 0.0)?; let k = e.tangent.take().ok_or_else(|| Failure::Fatal(Error::Internal("tangent".into())))?; let [kff] = self.problem.red.reduce([&k])?; let u = factor_raw(&kff, self.limits).map(|_| ()); @@ -670,7 +694,7 @@ impl Run<'_> { /// factored (the increment then starts from `q0`). fn tangent_predictor(&mut self, q0: &[f64], from: f64, to: f64) -> std::result::Result, f64)>, Failure> { let red = self.problem.red; - let mut e = self.problem.evaluate(&red.expand_scaled(q0, from), to, true, &self.committed, &self.multipliers)?; + let mut e = self.problem.evaluate(&red.expand_scaled(q0, from), to, true, &self.committed, &self.contact, &self.reference, self.band)?; let k = e.tangent.take().ok_or_else(|| Failure::Fatal(Error::Internal("tangent".into())))?; let rf: Vec = free_residual(red, &e).iter().zip(prescribed_force(red, &k, to - from)).map(|(a, b)| a - b).collect(); let mut dq = rf.clone(); @@ -756,6 +780,8 @@ impl Run<'_> { if sets.0 == Some(set) && set != sets.1 { flips += 1; if flips >= contact::CHATTER_LIMIT { + // From now on the nodes about to touch stay in the tangent. + self.band = contact::ENGAGED_BAND; return Err(Failure::Cutback("the contact status alternates between two states (chattering)".into())); } } @@ -816,9 +842,9 @@ impl Run<'_> { .problem .contacts .iter() - .zip(&self.multipliers) + .zip(&self.contact) .zip(&c.evaluation.contact) - .fold(0.0f64, |m, ((i, l), s)| m.max(i.augmentation_error(l, s))); + .fold(0.0f64, |m, ((i, l), s)| m.max(i.augmentation_error(&l.multipliers, s))); if error <= contact::GAP_TOLERANCE { return Ok(c); } @@ -832,7 +858,9 @@ impl Run<'_> { *count += 1; self.augmentations += 1; self.augmentation = *count; - self.multipliers = c.evaluation.contact.iter().map(|s| s.pressure.clone()).collect(); + for (l, s) in self.contact.iter_mut().zip(&c.evaluation.contact) { + l.multipliers.clone_from(&s.pressure); + } // From the converged state, measured against the increment's start and its first // correction like the increment itself. let start = std::mem::take(&mut c.q); @@ -1151,7 +1179,7 @@ fn solve_on(model: &Model, step: &str, limits: &Limits, notes: Vec) -> R .collect::>>()?; // The undeformed state first: an element that is bad before anything moves is an error in // the model, not a reason to cut increments back. - problem.evaluate(&vec![0.0; mesh.dofs()], 0.0, false, &committed, &[])?; + problem.evaluate(&vec![0.0; mesh.dofs()], 0.0, false, &committed, &[], &[], 0.0)?; let path_dependent = problem.materials.iter().any(MaterialModel::is_path_dependent); let mut run = Run { problem: &problem, @@ -1165,7 +1193,9 @@ fn solve_on(model: &Model, step: &str, limits: &Limits, notes: Vec) -> R committed, path_dependent, undeformed: None, - multipliers: problem.contacts.iter().map(|c| vec![0.0; c.nodes.len()]).collect(), + contact: problem.contacts.iter().map(Interface::committed).collect(), + reference: vec![0.0; mesh.dofs()], + band: 0.0, augmentations: 0, augmentation: 0, converged_scale: 0.0, @@ -1232,14 +1262,14 @@ fn solve_on(model: &Model, step: &str, limits: &Limits, notes: Vec) -> R Err(Failure::Fatal(e)) => return Err(e), } } - let multipliers = run.multipliers.clone(); + let saved = run.contact.clone(); let mut augmentations = 0; run.augmentation = 0; let solved = run .increment(&q, predictor.as_deref(), predicted_energy, target_lambda, change.abs(), &mut iterations) .and_then(|c| run.augment(c, &q, target_lambda, change.abs(), &mut iterations, &mut augmentations, &mut warnings)); if solved.is_err() { - run.multipliers = multipliers; + run.contact = saved; } match solved { Ok(c) => { @@ -1261,6 +1291,10 @@ fn solve_on(model: &Model, step: &str, limits: &Limits, notes: Vec) -> R q = c.q; keep_frame(&mut history, lambda, &red.expand_scaled(&q, lambda), &mut warnings); run.committed.clone_from(&c.evaluation.history); + for (l, s) in run.contact.iter_mut().zip(&c.evaluation.contact) { + l.commit(s); + } + run.reference = red.expand_scaled(&q, lambda); run.converged_scale = run.converged_scale.max(norm_inf(&c.evaluation.external)).max(norm_inf(&c.evaluation.internal)); last = Some(c.evaluation); fast = if count <= FAST_ITERATIONS { fast + 1 } else { 0 }; @@ -1300,7 +1334,7 @@ fn solve_on(model: &Model, step: &str, limits: &Limits, notes: Vec) -> R let u = red.expand_scaled(&q, lambda); let evaluation = match last { Some(e) => e, - None => problem.evaluate(&u, lambda, false, &run.committed, &run.multipliers)?, + None => problem.evaluate(&u, lambda, false, &run.committed, &run.contact, &run.reference, 0.0)?, }; if let Some((at, negative)) = history .increments diff --git a/crates/solver/tests/contact.rs b/crates/solver/tests/contact.rs index bc52305..9f93692 100644 --- a/crates/solver/tests/contact.rs +++ b/crates/solver/tests/contact.rs @@ -403,3 +403,63 @@ fn hertz_sphere_on_a_flat() { assert!((found - peak).abs() < 0.05 * peak, "p0 {found} vs {peak}"); assert!(before <= a && a <= next, "a {a} not in ({before}, {next})"); } + +/// Two blocks, the upper one pressed down by `down` and moved sideways by `side` at its top +/// (both in proportion over the step), with Coulomb friction `friction`; the reactions at the +/// top (normal, tangential) and the result. +fn sheared(down: f64, side: f64, friction: Friction, finite: bool) -> (f64, f64, StaticResult) { + let mut b = Blocks::default(); + let low = b.block(ElementKind::Hex8, [3, 3, 1], |p| p); + let up = b.block(ElementKind::Hex8, [2, 2, 1], |p| [p[0], p[1], 1.0 + p[2]]); + b.define(1000.0, 0.3, NonlinearStatic { nonlinear_geometry: finite, initial_increment: 0.25, ..NonlinearStatic::default() }); + b.face_set("lower_top", &low, |x| near(x[2], 1.0)); + b.face_set("upper_bottom", &up, |x| near(x[2], 1.0)); + let mut c = Contact::new(ContactSurfaces::Pair { main: "lower_top".into(), secondary: "upper_bottom".into() }); + c.friction = friction; + b.model.add_contact("interface", c).unwrap(); + b.bc("base", [Some(0.0), Some(0.0), Some(0.0)], |x| near(x[2], 0.0)); + b.bc("top", [Some(side), Some(0.0), Some(-down)], |x| near(x[2], 2.0)); + let r = solve_nonlinear(&b.model, Some("Step"), &Limits::default()).unwrap(); + report(&r); + let sum = |k: usize| -> f64 { b.nodes(|x| near(x[2], 2.0)).iter().map(|l| r.reaction[r.nodes.binary_search(l).unwrap()][k]).sum() }; + (-sum(2), sum(0), r) +} + +/// Coulomb's law across a frictional interface: a tangential demand below μ N sticks (the blocks +/// shear as one, nothing slips), above it the interface slips and carries exactly μ N; with a +/// kinetic coefficient the slipping interface carries μ_k N. +#[test] +fn coulomb_friction_sticks_below_the_cone_and_slips_on_it() { + for finite in [false, true] { + let mu = 0.3; + // Stuck as one body, the shear the top displacement asks for: a frictionless trial gives + // the normal force, a stuck one with a very large μ the tangential force. + let (n, t, _) = sheared(0.01, 0.001, Friction::coulomb(5.0), finite); + let ratio = t / n; + eprintln!("finite {finite}: stuck N {n:.4} T {t:.4} ratio {ratio:.4}"); + // Well below the cone (the shear concentrates at the edges of the interface, which slip + // first): sticks, the same forces as glued. + let below = 0.4 * mu / ratio; + let (n1, t1, r1) = sheared(0.01, 0.001 * below, Friction::coulomb(mu), finite); + assert!((t1 / n1 - 0.4 * mu).abs() < 1e-3 * mu, "finite {finite}: {} vs {}", t1 / n1, 0.4 * mu); + // The middle of the interface sticks (its corners, where the shear concentrates, may slip). + // (The upper block's 3 × 3 bottom nodes come first, in x then y order: the fifth.) + assert_eq!(r1.contacts[0].status[4], Status::Stick, "finite {finite}: {:?}", r1.contacts[0].status); + // Above it: slips, carrying μ N. + let (n2, t2, r2) = sheared(0.01, 0.001 * 2.0 * mu / ratio, Friction::coulomb(mu), finite); + // With finite sliding the upper block turns a little as it slides, and the contact acts + // along its face's normal, a few thousandths of a radian from the vertical. + let tol = if finite { 0.03 } else { 1e-3 }; + assert!((t2 - mu * n2).abs() < tol * mu * n2, "finite {finite}: T {t2} vs μN {}", mu * n2); + assert!(r2.contacts[0].status.iter().all(|s| *s == Status::Slip), "finite {finite}: {:?}", r2.contacts[0].status); + let shear: f64 = r2.contacts[0].shear.iter().map(|s| s[0]).zip(&r2.contacts[0].pressure).map(|(s, p)| s / p).sum::() + / r2.contacts[0].nodes.len() as f64; + // Node by node the nodal forces' ratio scatters a little about μ (the smoothing of the + // nodal pressure and shear differs at the interface's edges); the mean is μ. + assert!((shear + mu).abs() < 0.05 * mu, "finite {finite}: shear over pressure {shear}"); + // Kinetic friction once slipping. + let decaying = Friction { static_coefficient: mu, kinetic: Some(0.2), decay: 1.0 }; + let (n3, t3, _) = sheared(0.01, 0.001 * 2.0 * mu / ratio, decaying, finite); + assert!((t3 - 0.2 * n3).abs() < tol * 0.2 * n3, "finite {finite}: T {t3} vs μk N {}", 0.2 * n3); + } +}