Kuwabara-Kono numerical dissipation: a new method to simulate granular matter

Abstract

A new method is introduced for the simulation of multiple impacts in granular media using the Kuwabara-Kono (KK) contact model, a nonsmooth (not Lipschitz continuous) extension of Hertz contact that accounts for viscoelastic damping. We use the technique of modified equations to construct time-discretizations of the nondissipative Hertz law matching numerical dissipation with KK dissipation at different consistency orders. This allows us to simulate dissipative impacts with good accuracy without including the nonsmooth KK viscoelastic component in the contact force. This tailored numerical dissipation is developed in a general framework, for Newtonian dynamical systems subject to dissipative forces proportional to the time-derivative of conservative forces. Numerical tests are performed for the simulation of impacts in Newton’s cradle and on alignments of alternating large and small balls. Resulting wave phenomena (oscillator synchronization, propagation of dissipative solitary waves, oscillatory tails) are accurately captured by implicit schemes with tailored numerical dissipation, even for relatively large time steps.

Nonlinear phenomena of contact in multibody systems dynamics: a review

“…It is assumed that materials are isotropic and that there are no plastic deformations [170]. This model was used by James and his co-authors to develop a method for simulating multiple impacts in granular mediums [171].…”

Traveling fronts in dissipative granular chains and nonlinear lattices

“…We also apply the scaling properties to the normalization of front solutions. Section 3.2 recalls dynamical variables introduced in [30] for the numerical simulation of the models, which are adapted to the limited smoothness of the viscoelastic contact force. In section 3.3, we simulate the KK model (case α = β = 3/2) with a constant pressure applied at one end of the chain, and show that underdamped or overdamped fronts are generated depending on parameter values.…”

“…In this range of parameters, the model presents a lack of smoothness because the map defining the right-hand side is discontinuous accross the hyperplanes x n−1 = x n for β = 1, and not Lipschitz continuous for β ∈ (1, 2), in particular for the KK model. This lack of regularity can induce a rather severe decrease of the order of convergence of some classical numerical time-integration schemes, as shown in [30], and may as well deteriorate the computation of wave profiles based on Newton-type methods. However, these difficulties can be circumvented to some extent by an appropriate choice of dynamical variables [30].…”

“…This lack of regularity can induce a rather severe decrease of the order of convergence of some classical numerical time-integration schemes, as shown in [30], and may as well deteriorate the computation of wave profiles based on Newton-type methods. However, these difficulties can be circumvented to some extent by an appropriate choice of dynamical variables [30]. Defining a generalized velocity w n througḣ…”

“…The KK model has the advantage of being derived from continuum mechanics (see [6,18] and references therein) and γ 0 is determined by the elastic and viscous coefficients of the material. Moreover, it was shown in [30] that the KK model approximates the numerical dissipation of a class of first order Runge-Kutta schemes (including the θ method) used to integrate nondissipative granular chains. Consequently, the KK model can be also useful in the context of nondissipative Hertzian chains, in order to evaluate the effect of numerical damping on the dynamics of nonstationary shock waves (a problem suggested in [39]).…”