Hyperbolic dynamic friction models for viscoelastic sliding and rolling contact

1 Introduction

• When a sliding substrate is compliant or both bodies deform, material advection undermines the premise that sliding friction can be represented by finite-dimensional ODE models.Existing bristle-based sliding models commonly use nonlinear ODEs for internal states, while rolling models are typically PDEs.Source: §1 · p. 2
• The work extends the purely elastic analysis in to viscoelastic contact pairs, where constitutive behavior enters the kinematic relationships for dissimilar bodies.The extension addresses both sliding and rolling friction for viscoelastic pairs.Source: §1 · p. 2
• The framework combines linear viscoelastic rheologies, including GKV and GM equivalents, with nonlinear bristle-state dynamics to derive semilinear PDEs for contact forces and internal states.The stated scope is classical derivative-based linear viscoelasticity and planar, localized contact configurations.Source: §1 · p. 2
• Sliding and rolling can share the same linear or semilinear hyperbolic PDE structure, making their traditional ODE-versus-PDE distinction an artificial modeling split.The framework ties equation form to contact kinematics and constitutive behavior rather than assigning a formulation by contact type.Source: §1 · pp. 1, 3
• The framework is not a complete predictive contact solver and excludes normal–tangential coupling, thermoviscoelastic effects, and nonlinear rheological behavior.These mechanisms are deliberately neglected to keep the analytical theory tractable.Source: §1 · p. 3

2.1 General structure of dynamic friction models

• The general model tracks relative bristle deformation, friction force, and internal states, with rheological and force-evolution functions specifying their dynamics.With advection, the state derivatives become material derivatives and can yield hyperbolic PDEs; different transport velocities complicate this extension.Source: §2.1 · pp. 4–5

2.2 Review of dynamic friction models

• The review restricts attention to sufficiently smooth dynamic friction laws, including Dahl, LuGre, FrBD, and Bouc-Wen formulations.The specific structure of the force-evolution function does not determine whether the dynamics are lumped or distributed.Source: §2.2 · p. 5, §2.1 · p. 5

2.2.1 Dahl and LuGre models

• Dahl and LuGre specify the force-evolution function using a positive-definite friction-coefficient matrix, micro-stiffness coefficients, and a regularized Euclidean norm.The norm regularization is continuous and often converges uniformly to the Euclidean norm as its regularization parameter tends to zero.Source: §2.2.1 · pp. 5–6

2.2.2 FrBD models

• FrBD replaces the elastic micro-stiffness term in the Dahl/LuGre form with the generated force, with implications for model passivity.A related matrix term may be replaced by a positive-semidefinite matrix when the steady-state friction force has the specified form.Source: §2.2.2 · p. 6

2.2.3 Bouc-Wen model

• The two-dimensional Bouc-Wen model specifies the force-evolution function using direction-specific terms and parameters β1, β2, and γ.Absolute values in the model may be regularized using a similar approach to that used for Dahl and LuGre.Source: §2.2.3 · p. 6

2.2.4 Frictionless contact

• The frictionless-contact function can be obtained by linearizing the Dahl and FrBD formulations near zero deformation, force, and relative velocity.The resulting expression agrees with linear regimes of other friction models.Source: §2.2.4 · p. 6

3 Viscoelastic contact pairs and viscoelasto-kinematics

• Linear viscoelastic constitutive equations for bristles provide the dynamics of relative bristle deflection and frictional force.The section extends from individual bristle elements to viscoelastic contact pairs.Source: §3 · p. 6, §3.1 · p. 7

3.1.1 Linear viscoelastic rheological models

• Generalised Kelvin-Voigt and Generalised Maxwell rheologies combine elastic and viscous elements and admit stable, dissipative state-space representations.The GM model places an elastic element in parallel with dissipative branches; the models use force input and deformation output.Source: §3.1.1 · p. 7

3.1.2 Viscoelastic contact pairs

• The contact-pair formulation allows the two bodies to have distinct stable linear bristle rheologies, while assuming coaxial elastic stiffness matrices in the main analysis.The coaxiality condition is automatically satisfied for isotropic materials or anisotropic materials sharing principal elastic directions, and is removed in Appendix A.3 without altering the main conclusions.Source: §3.1.2 · p. 8

3.2.1 Apparent contact area and Eulerian derivatives

• Contact areas are treated in a contact-fixed frame, and Eulerian material derivatives account for quantities advected with each body's transport velocity.The apparent contact area is assumed compact and regular closed; the two bodies' transport velocities may differ.Source: §3.2.1 · pp. 9–10

3.2.2 Viscoelasto-kinematic equations

• The viscoelasto-kinematic derivation combines separate body coordinates using a common interface coordinate and assumes small displacements and displacement gradients.The first-order approximation neglects products of O(ϵ) Jacobian corrections with bounded spatial gradients and does not resolve layers whose gradients grow as O(ϵ^-1).Source: §3.2.2 · pp. 10–11
• Equations (25), (26), and (29) govern sliding and rolling contact, and the bristle force is the natural state variable in the viscoelastic setting.Writing the dynamics in terms of bristle displacement instead would introduce higher-order coupling terms with internal variables.Source: §3.2.2 · pp. 11–12
• The shared coordinate is a stiffness-weighted average of the two body coordinates plus a viscoelastic term, and the paper identifies this definition as novel.In the elastic case, the viscoelastic term vanishes and the coordinate reduces to the definition in.Source: §3.2.2 · p. 12

4.1 Model equations

• The governing equations are divided into PDEs for bristle-force and internal-variable evolution and equilibrium equations for computing global forces and moments.The evolution equations depend on rigid relative velocity input.Source: §4.1 · p. 12

4.1.1 Governing PDEs

• The proposed equations are not simple transport equations, yet remain hyperbolic and derive their advection terms from material transport through the contact region.When transport disappears, the distributed nature of the friction dynamics can be reduced.Source: §4.1.1 · p. 13

4.1.2 Forces and moments (global equilibrium)

• Global forces and moments follow from the PDE solution, and the formulation makes distributed or lumped dynamics emerge from material transport rather than a modelling assumption.The total vertical-moment formula in Eq. (38) is described as novel, whereas Eq. (39) is standard.Source: §4.1.2 · pp. 13–14

4.2 Hyperbolicity and Riemann variables

• The PDE system retains rolling contact’s hyperbolic character, and a change to Riemann variables restores a transport-like form.Source: §4.2 · p. 14

4.2.1 Hyperbolicity

• Under Assumptions 3.1 and 3.2, the PDE system is symmetric hyperbolic; the stated proof setup assumes a time-independent contact domain.The authors conclude that viscoelasticity does not alter the hyperbolic nature of the frictional processes.Source: §4.2.1 · p. 14

4.2.2 Riemann variables

• In Riemann variables, force components propagate along characteristic directions set by the contact pair’s principal elastic interaction modes.The transformed system can be interpreted as four vector-valued transport equations with different transport velocities.Source: §4.2.2 · p. 15

4.3 Boundary and initial conditions

• Boundary conditions enforce continuity of principal stresses at the contact transition and ensure dissipativity; initial conditions can be prescribed in Riemann variables.The system’s well-posedness is stated for time-invariant contact areas and velocities.Source: §4.2.2 · p. 15, §4.3 · p. 16

5.1 Sliding friction

• Sliding couples hyperbolic PDEs for advected substrate bristles with spatially local ODEs for non-advected block bristles, and the constant-speed system is well-posed.The local ODE dynamics vary with position, unlike spatially uniform finite-dimensional friction models.Source: §5.1 · pp. 17–18

5.1.1 Steady sliding

• As s decreases, steady sliding forces increase; adding substrate dissipative branches reduces interface force through spatial relaxation, while longer relaxation times increase force.The steady response is independent of block rheology.Source: §5.1.1 · p. 19, Figure 5 · p. 19

5.1.2 Transient sliding

• At the leading edge during transient sliding, block-bristle deformation converges exponentially to zero unless its internal states start at equilibrium; lower s raises force and longer relaxation times speed convergence.Block rheology contributes during transients, unlike in the steady response.Source: §5.1.2 · pp. 19, 21

5.2 Rolling friction

• Rolling retains full PDE advection for all positive rolling and forward velocities, and the constant-velocity system has existence and uniqueness properties.Unlike sliding, rolling does not lose advection in the modeled equations.Source: §5.2 · pp. 21, 23

5.2.1 Steady rolling

• In steady rolling, longitudinal force depends on relative velocity, rolling speed, and s, vanishes when relative velocity is zero, and shows only minor influence from higher-order rheologies.The force depends on the kinematic difference between rolling and translation because both bodies advect through the contact patch.Source: §5.2.1 · p. 23

5.2.2 Transient rolling

• After a rolling relative-velocity step, decreasing s from 0.8 to 0.4 produces slightly higher longitudinal forces over shorter time scales.The example uses Vr = 16 m s−1, Vx = 12.8 m s−1, and vx = −3.2 m s−1.Source: §5.2.2 · pp. 23–24

• The framework combines linear viscoelastic rheologies with bristle-based dynamic friction laws to derive semilinear PDEs for friction forces, bristle deformation, and interface state variables.Source: §6 · pp. 24–25
• Dynamic friction models extend from lumped ODEs to distributed equations that account for material advection, including contact pairs with differing elastic and viscoelastic properties.Source: §6 · pp. 24–25
• Under small deformations and displacement gradients, the governing equations retain a symmetric or symmetrisable hyperbolic structure and can be recast as transport equations along elastic-interaction characteristics.Source: §6 · p. 25
• In sliding, steady response is largely governed by elastic deformation partition, while higher-order viscoelastic rheologies mainly affect transient behaviour and dissipation.Source: §6 · p. 25
• In rolling, advection in both bodies preserves distributed dynamics and produces force-generation mechanisms distinct from sliding, despite mathematical analogies between the regimes.Source: §6 · p. 25
• The framework’s future development includes nonlinear viscoelasticity, viscoplasticity, efficient large-scale numerical schemes, experimental validation, and contact configurations with spin.Source: §6 · p. 26

• The FASTEST project (Reg. no. 2023-06511), funded by the Swedish Research Council, financially supported this research.Source: Funding declaration · p. 26

• L.R. is the sole author and contributor to the manuscript.Source: Author Contribution declaration · p. 26

• The appendix shows hyperbolicity is retained for variable contact areas, proves well-posedness for time-independent areas and velocities, and discusses extension to the non-coaxial case.Source: A Technical results · p. 30

A.1 Hyperbolicity in the case of a time-varying contact area

• A level-preserving smooth diffeomorphism from the time-varying domain to a cylinder preserves the PDE system’s hyperbolic character.The transformed system is symmetric hyperbolic under Assumptions 3.1 and 3.2.Source: §A.1 · p. 31
• The transformed boundary conditions remain valid on a time-varying contact area, and the hyperbolic system can be reformulated in Riemann variables.The Riemann-variable reformulation and boundary-condition result are stated by analogy with the time-invariant case.Source: §A.1 · p. 31

A.2 Well-posedness

• The transformed system admits a unique mild solution when its nonlinear term is continuous in time and uniformly Lipschitz, for any L2 initial condition.With a continuously differentiable nonlinear term and initial data in D(A), the solution is classical.Source: §A.2 · p. 32
• For LuGre and FrBD models, global solutions are established in restricted settings, including cases without spin or with approximated relative velocity.The stated cases include non-regularised LuGre and regularised FrBD friction laws; smoother laws may allow more regular solutions under additional conditions.Source: §A.2 · p. 32

A.3 Extension to non-coaxial elastic stiffness matrices

• Without coaxiality, the reformulated PDE system remains symmetrisable hyperbolic but generally requires modified equations and lacks the orthogonal characteristic representation.A common positive-definite symmetriser exists even though the transport matrices are not generally symmetric.Source: §A.3 · pp. 32–33
• In the non-coaxial case, S1 and S2 still partition deformation, while their transposes determine force transport and can be simultaneously diagonalised.The transposed matrices commute and are self-adjoint in the K0^-1-weighted inner product.Source: §A.3 · pp. 33–34

B FrBDn+1-GKV model for viscoelastic contact pairs

• The FrBDn+1-GKV model is presented for viscoelastic contact pairs, with derivation details deferred and passivity treated separately.Its equations assume time-independent contact areas and velocities, with extension to time-varying cases described as immediate.Source: B FrBDn+1-GKV model for viscoelastic… · p. 34, §B.1 · p. 34

B.1 Model equations

• The model equations specify internal variables and normalized stiffness, damping, and relaxation-time matrices for the two bodies.Under stated coefficient and input-velocity regularity conditions, the equations can be recast to establish mild or classical well-posedness.Source: §B.1 · p. 35

B.2 Passivity

• Under Assumption 3.1 and the stated pressure inequalities, the FrBDn+1-GKV contact model is passive with an energy-based storage function.The result extends to the non-coaxial case under analogous assumptions and appropriate boundary conditions.Source: §B.2 · pp. 36–37
• When transport velocities are constant in space, a constant pressure distribution satisfies the passivity inequalities; decreasing pressure along rolling can make storage a Lyapunov function.The pressure conditions may be verified simultaneously when the two transport velocities have the same direction.Source: §B.2 · p. 37

C.1 Semi-discrete formulation

• The spatially discretized equations are integrated in MATLAB with adaptive ode15s, while the maximum time step is tied to cell size and the largest transport velocity.The stated tolerances are RelTol = 10^-9 and AbsTol = 10^-11; the step restriction resolves material transport across a spatial cell.Source: §C.1 · p. 38

C.2 Numerical convergence study

• The convergence study examines both steady-state and transient problems using the model parameters from Tables 1 and 2.Source: §C.2 · p. 38

C.2.1 Verification in the purely elastic limit for the stationary solutions

• The verification uses relative spatial L2 and L∞ errors and reports q as the observed convergence order in the relative L2 norm.Source: Table 3 · p. 39
• Numerical stationary bristle forces are compared with exact elastic solutions for sliding and rolling, with faster convergence reported for rolling.The sliding results are already satisfactory at N = 101.Source: §C.2.1 · p. 39, Figure 10 · p. 39
• Rolling’s observed order qk rises from 0.985 to 0.996, while sliding’s rises from 0.772 to 0.936, indicating convergence toward first order.The trends are consistent with the first-order upwind spatial discretisation.Source: §C.2.1 · p. 39

C.2.2 Transient convergence analysis

• Because no exact full viscoelastic transient solution is available, the study uses coupled space-time refinement, reducing the maximum integration step proportionally to the spatial step.Force histories are sampled at common output times and compared with the Nref = 401 solution.Source: §C.2.2 · pp. 39–40
• Figure 11 tracks elastic steady-state and transient convergence errors across N = 51, 101, 201, and 401 at s = 0.4.The finest-grid transient difference is zero by construction and is omitted from the logarithmic plot.Source: Figure 11 · p. 40, §C.2.2 · p. 41
• The refinement analysis reports relative L2- and L∞-time errors and relative differences in peak and terminal resultant forces.Table 4 lists these quantities for s = 0.4; the N = 401 self-comparison entries are zero by definition.Source: Table 4 · p. 40, §C.2.2 · p. 40
• The numerical longitudinal forces converge rapidly under refinement for both sliding and rolling contact.The relative L2-time differences also decrease systematically as the grid is refined.Source: §C.2.2 · pp. 40–41