Hertz Model

The Hertz contact model consists of a nonlinear formulation based on an approximation of the theory of Mindlin and Deresiewicz [Mindlin1953]. This model can be installed at both ball-ball and ball-facet contacts, and is referred to in commands and FISH by the name hertz.

Introduction

The Hertz contact model produces both normal and shear forces based on the theoretical analysis of the deformation of smooth, elastic spheres in frictional contact. Viscous dashpots may be added for further energy dissipation in simulations involving impact. Both components transmit only a force.

Behavior Summary

The Hertz contact model consists of a nonlinear Hertzian component based on an approximation of the theory of Mindlin and Deresiewicz [Mindlin1953]. As discussed in [Cundall1988a], the current implementation does not reproduce the continuous nonlinearity in shear. Instead, the initial shear modulus is used, but it depends on normal force. This formulation presumes a finite-sized contact area that is small relative to the radii of curvature of the contacting pieces. The interface does not resist relative rotation, meaning that the contact moment equals zero (\(\mathbf{M_{c}} \equiv 0\)). The contact force is resolved into Hertzian and dashpot components (\(\mathbf{F_{c}} =\mathbf{F^{h}} +\mathbf{F^{d}}\)). A surface gap, \(g_{s}\), is defined as the difference between the contact gap \(g_{c}\) and the reference gap \(g_{r}\) so that when the reference gap is zero, the notional surfaces coincide with the piece surfaces. The contact is active if and only if the surface gap is less than or equal to zero; the force-displacement law is skipped for inactive contacts.

common\contactmodel\hertz\doc\manual\cmhertz_fig1.png

Figure 1: Behavior and rheological components of the Hertz model.

Note

The Hertz model implementation has changed from PFC 4.0 as follows:

  1. Ball-Facet Contact Treatment: In PFC 4.0, the wall in a ball-wall contact was conceptually replaced by an identical ball positioned symmetrically about the wall location so that the value of the overlap used in the force-displacement law was set equal to twice the physical ball-wall overlap. In PFC 5.0 and PFC 6.0, the force is resolved for a ball-facet contact using the facet curvature of 0 (or, alternatively, the facet radius of \(\infty\)). This means that the normal force at a ball-facet contact is half the value calculated in PFC 4.0.
  2. 2D Formulation: The theoretical background of the Hertz model in PFC is only applicable to the case of 3D ball-ball and ball-facet contacts. Users may however use the Hertz contact model in PFC2D and select the value of the power-law exponent \(\alpha_h\) (see below).
  3. The power-law exponent \(\alpha_h\) (see description below) can be modified from its default value of 3/2. For instance, [Hunt1975] suggest a value in the range [1.0,1.5] for contacts between coaxial cylinders.
  4. A new mode (\(M_s\)) allows for the shear force to be diminished upon normal unloading to prevent spurious energy creation (see [Ellata1996] or [Agnolin2007a] for details).

Note

Revision 28 included the following additional modifications:

  1. The property inheritance mechanism that derives the contact model elastic parameters from the parameters of the contacting pieces has been changed to conform to [Mindlin1953], following the equations below . Before this modifications, the contact effective shear modulus and Poisson ratio were taken as the average of the contacting pieces shear modulus and Poisson ratio, respectively. Note that this change may change the output for models that used the property inheritance mechanism with the hertz model. Please refer to the section “Usage and Verification Examples” below for an example demonstrating how to override the property inheritance mechanism to reproduce the previous behavior.
  2. A new property, \(\alpha_d\), has been added to allow the dashpot force to be proportional to the current surface gap raised to the power \(\alpha_d\). This parameter defaults to 0.0 for backwards compatibility, but can be set to a non-zero value. Following the dimensional analysis of [Ramírez1999], \(\alpha_d = (\alpha_h-1) / 2\) should be selected to obtain a constant restitution coefficient at low impact velocities. The current model with \(\alpha_h = 3/2\) and \(\alpha_d = 1/4\) satisfies this condition and corresponds to that of [Tsuji1992].
  3. A new property, \(g_r\), has been added to allow the user to specify a non-zero reference gap.

Activity-Deletion Criteria

A contact with the Hertz model is active if and only if the surface gap (\(g_{s} =g_{c} -g_{r}\)) is less than or equal to 0. As discussed in “Contact Resolution,” force-displacement calculations are skipped for inactive contacts.

Force-Displacement Law

The force-displacement law for the Hertz model updates the contact force and moment:

(1)\[\mathbf{F_{c}} =\mathbf{F^{h}} +\mathbf{F^{d}} ,\quad \mathbf{M_{c}} \equiv 0\]

where \(\mathbf{F^{h}}\) is the nonlinear Hertz force and \(\mathbf{F^{d}}\) is the dashpot force. Both forces are resolved into normal and shear forces:

(2)\[\mathbf{F^{h}} =-F_{n}^{h} {\kern 1pt} \hat{\mathbf{n}}_\mathbf{c} +\mathbf{F_{s}^{h}} ,\quad \mathbf{F^{d}} =-F_{n}^{d} {\kern 1pt} \hat{\mathbf{n}}_\mathbf{c} +\mathbf{F_{s}^{d}}.\]

where \(F_{n}^{\left\{h,d\right\}} >0\) is tension. Both shear forces lie on the contact plane, and the Hertz shear force is expressed in the contact plane coordinate system:

(3)\[\mathbf{F_{s}^{h}} =F_{ss}^{h} {\kern 1pt} \hat{\mathbf{s}}_\mathbf{c} +F_{st}^{h} {\kern 1pt} \hat{\mathbf{t}}_\mathbf{c} \; \; \left({\rm 2D\; model:\; }F_{ss}^{h} \equiv 0\right).\]

The force-displacement law for the Hertz model consists of the following steps:

  1. Update the Hertz normal force:

    (4)\[\begin{split}F^{h}_n = \begin{cases} - h_n {|g_c|}^{\alpha_h} & \mbox{,} g_c < 0.0 \\ 0.0 & \mbox{, otherwise} \end{cases}\end{split}\]

    The coefficient \(h_n\) depends on the geometrical and mechanical properties of the contacting pieces as

    (5)\[h_n = \frac{2 G \sqrt{2\bar{R}}}{3(1-\nu)}\]

    where \(G\), \(\nu\), and \(\bar{R}\) are the effective shear modulus, Poisson’s ratio, and contact radius, respectively. The effective shear modulus and Poisson’s ratio can be input directly, or inherited from the contacting pieces as detailled below.

    The effective radius of the contact is computed via the radii of the contacting piece surfaces (where the radius of a facet is \(\infty\)) via

    (6)\[\frac{1}{\bar{R}} = \frac{1}{2}\left(\frac{1}{R^{(1)}} + \frac{1}{R^{(2)}}\right)\]

  2. Update the Hertz model shear force as follows. Compute a trial shear force:

    (7)\[\mathbf{F^{*}_s} = \mathbf{({F}^{h}_s)_0} + k_s \Delta \pmb{θ}_\mathbf{s}\]

    where \(\mathbf{({F}^{h}_s)_0}\) is the Hertz shear force at the beginning of the timestep and \(\Delta \pmb{θ}_\mathbf{s}\) is the relative shear increment in this equation of the “Contact Resolution” section. In the expression above, \(k_s\) is the initial tangent shear stiffness, which depends on the current normal force as

    (8)\[k_s = \frac{2(1-\nu)}{2-\nu} \alpha_h \left(h_n\right)^{1/\alpha_h} \left(F^{h}_n\right)^{(\alpha_h-1)/\alpha_h}\]

    The value of \(\mathbf{({F}^{h}_s)_0}\) depends on the scaling mode (\(M_s\)):

    (9)\[\begin{split}\mathbf{({F}^{h}_s)_0} = \begin{cases} \frac{k_s}{(k_s)_0} \mathbf{(F^{h}_s)_0} & \mbox{, if } M_s = 1 \mbox{ and } F^{h}_n < (F^{h}_n)_0 \\ \\ \mathbf{(F^{h}_s)_0} & \mbox{, otherwise} \end{cases}\end{split}\]

    where \((F^{h}_n)_0\), \(\mathbf{(F^{h}_s)_0}\), and \((k_s)_0\) are the normal force, the shear force, and the tangent shear stiffness computed at the beginning of the timestep, respectively. Compute the shear strength:

    (10)\[F^{\mu}_s = \mu F^{h}_n.\]

    Update the Hertz shear force:

    (11)\[\begin{split}\mathbf{F_{s}^{h}} = \begin{cases} \mathbf{F_{s}^{*}} & \mbox {, } \left\| \mathbf{F_{s}^{*}} \right\| \le F_{s}^{\mu} \\ F_{s}^{\mu} \left(\mathbf{F_{s}^{*}} / \left\| F_{s}^{*} \right\| \right) & \mbox {, otherwise.} \end{cases}\end{split}\]

    Update the slip state:

    (12)\[\begin{split}s= \begin{cases} {\rm true} &\mbox{, } \left\| \mathbf{F_{s}^{h}} \right\| =F_{s}^{\mu } \\ {\rm false} & \mbox{, otherwise.} \end{cases}\end{split}\]

    If the slip state is true, then the contact is sliding. Whenever the slip states changes, the slip_change callback event occurs.

  3. Update the dashpot normal force based on the dashpot behavior mode:

    (13)\[\begin{split}\begin{array}{l} {F_{n}^{d} = \left\{\begin{array}{l} {\qquad F^{*} ,\qquad M_{d} =\left\{0,2\right\}{\rm \; (full\; normal)}} \\ {\qquad \min \left(F^{*} ,-F_{n}^{l} \right), \qquad M_{d} =\left\{1,3\right\}{\rm \; (no-tension\; normal)}} \end{array} \right. } \\ {{\rm with}\qquad F^{*} =\left(2\beta _{n} \sqrt{m_{c} k_{n} } \right) {|g_c|}^{\alpha_d} \dot{\delta }_{n} } \\ {\qquad k_{n} = {\qquad \alpha_h h_n \delta_n^{\alpha_h-1}}} \\ {\qquad m_{c} =\left\{\begin{array}{l} {\qquad \frac{m^{\left(1\right)} m^{\left(2\right)} }{m^{\left(1\right)} +m^{\left(2\right)} } ,\qquad {\rm ball-ball}} \\ {\qquad m^{\left(1\right)} {\rm ,}\qquad {\rm ball-facet}} \end{array}\right. } \end{array}\end{split}\]

    where \(k_{n}\) is the normal tangent stiffness, \(m^{\left(b\right)}\) is the mass of body (\(b\)), \(\alpha_d\) is a non-dimensionnal exponent, and \(\dot{\delta }_{n}\) is the relative normal translational velocity of this equation of the “Contact Resolution” section. For the dashpot normal-behavior mode of no-tension, \(F_{n}^{d}\) is capped to ensure that the total normal force (\(F_{n}^{l} +F_{n}^{d}\)) does not become tensile.

  4. Update the dashpot shear force based on the dashpot behavior mode:

    (14)\[\begin{split}\begin{array}{l} {\mathbf{F_{s}^{d}} =\left\{\begin{array}{l} {\qquad \left(2\beta _{s} \sqrt{m_{c} k_{s} } \right) {|g_c|}^{\alpha_d} \dot{\pmb{δ}}_\mathbf{s} ,\qquad s={\rm false\; \; or}\qquad M_{d} =\left\{0,1\right\}{\rm \; (full\; shear)}} \\ {\qquad 0,\qquad s={\rm true\; \; and}\qquad M_{d} =\left\{2,3\right\}{\rm \; (slip-cut)}} \end{array}\right. } \\ {{\rm with}\qquad m_{c} =\left\{\begin{array}{l} {\qquad \frac{m^{\left(1\right)} m^{\left(2\right)} }{m^{\left(1\right)} +m^{\left(2\right)} } ,\qquad {\rm ball-ball}} \\ {\qquad m^{\left(1\right)} {\rm ,}\qquad {\rm ball-facet}} \end{array}\right. } \end{array}\end{split}\]

    where \(\dot{\pmb{δ}}_\mathbf{s}\) is the relative shear translational velocity of this equation of the “Contact Resolution” section. For the dashpot shear-behavior mode of slip-cut, \(\mathbf{F_{s}^{d}}\) is set to zero if the linear spring is sliding (\(s\) = true).

Energy Partitions

The Hertz model defines three energy partitions:

  • strain energy \(E_{k}\), stored in the springs;
  • slip work \(E_{\mu}\), defined as total energy dissipated by frictional sliding; and
  • damp work \(E_{\beta}\), defined as the total energy dissipated by the dashpots.

If energy tracking is activated, these energy partitions are updated as follows.

  1. Update the strain energy:

    (15)\[E_{k} = \frac{\alpha_h}{(\alpha_h+1)} \left(\frac{\left(F^{h}_n\right)^2}{k_n}\right) + \frac{1}{2} \frac{\|\mathbf{F^{h}_s)}\|^2}{k_s}\]

  2. Update the slip energy:

    (16)\[\begin{split}\begin{array}{l} {E_{\mu } :=E_{\mu } - {\tfrac{1}{2}} \left(\left(\mathbf{F_{s}^{h}} \right)_{o} +\mathbf{F_{s}^{h}} \right)\cdot \Delta \pmb{δ}_\mathbf{s}^{\mu } } \\ {{\rm with}\qquad \Delta \pmb{δ}_\mathbf{s}^{\mu } =\Delta \pmb{δ}_\mathbf{s} -\Delta \pmb{δ}_\mathbf{s}^{k} =\Delta \pmb{δ}_\mathbf{s} -\left(\frac{\mathbf{F_{s}^{h}} -\left(\mathbf{F_{s}^{h}} \right)_{o} }{k_{s} } \right)} \end{array}\end{split}\]

    where \(\left(\mathbf{F_{s}^{h}} \right)_{o}\) is the shear force at the beginning of the timestep, and the relative shear-displacement increment of this equation of the “Contact Resolution” section has been decomposed into an elastic \(\left(\Delta \pmb{δ}_\mathbf{s}^{k} \right)\) and a slip \(\left(\Delta \pmb{δ}_\mathbf{s}^{\mu } \right)\) component. The dot product of the slip component and the average shear force occurring during the timestep gives the increment of slip energy.

  3. Update the dashpot energy:

    (17)\[E_{\beta } :=E_{\beta } - \mathbf{F^{d}} \cdot \left(\dot{\pmb{δ} }{\kern 1pt} {\kern 1pt} \Delta t\right)\]

    where \(\dot{\pmb{δ} }\) is the relative translational velocity of this equation of the “Contact Resolution” section.

Additional information, including the keywords by which these partitions are referred to in commands and FISH, is provided in the table below.


Table 1: Hertz Model Energies
Keyword Symbol Description Range Accumulated
Hertz Group:
energy‑strain \(E_{k}\) strain energy \([0.0,+\infty)\) NO
energy‑slip \(E_{\mu}\) total energy dissipated by slip \([0.0,+\infty)\) YES
Dashpot Group:
energy‑dashpot \(E_{\beta}\) total energy dissipated by dashpots \([0.0,+\infty)\) YES

Properties

The properties defined by the Hertz contact model are listed in the table below for a concise reference; see the “Contact Properties” section for a description of the information in the table columns. The mapping from the surface inheritable properties to the contact model properties is also discussed below.

Table 2: Hertz Model Properties
Keyword Symbol Description Type Range Default Modifiable Inheritable
hertz Model name
Hertz Group:
rgap \(g_r\) Reference gap [length] FLT \(\mathbb{R}\) 0.0 YES NO
hz_shear \(G\) Shear modulus [stress] FLT \([0.0,+\infty)\) 0.0 YES YES
hz_poiss \(\nu\) Poisson’s ratio [-] FLT \((-1.0,0.5]\) 0.0 YES YES
fric \(\mu\) Friction coefficient [-] FLT \([0.0,+\infty)\) 0.0 YES YES
hz_mode \(M_s\) Shear-force scaling mode [-] INT \(\{0;1\}\) 0 YES NO
    \(\;\;\;\;\;\;\begin{cases} \mbox{0: no scaling} \\ \mbox{1: scaling is active} \end{cases}\)          
hz_alpha \(\alpha_h\) Exponent [-] FLT \([0.0,+\infty)\) 1.5 YES NO
hz_slip \(s\) Slip state [-] BOOL {false,true} false NO N/A
hz_force \(\mathbf{F^{h}}\) Hertz force (contact plane coord. system) VEC \(\mathbb{R}^3\) \(\mathbf{0}\) YES NO
    \(\left( -F_n^h,F_{ss}^h,F_{st}^h \right) \quad \left(\mbox{2D model: } F_{ss}^h \equiv 0 \right)\)          
Dashpot Group:
dp_nratio \(\beta_n\) Normal critical damping ratio [-] FLT \([0.0,1.0]\) 0.0 YES NO
dp_sratio \(\beta_s\) Shear critical damping ratio [-] FLT \([0.0,1.0]\) 0.0 YES NO
dp_alpha \(\alpha_d\) Exponent [-] FLT \([0.0,+\infty)\) 0.0 YES NO
dp_mode \(M_d\) Dashpot mode [-] INT {0;1;2;3} 0 YES NO
    \(\;\;\;\;\;\;\begin{cases} \mbox{0: always active} \\ \mbox{1: no tension} \\ \mbox{2: no shear damping if s=1} \\ \mbox{3: combine modes 1 \& 2} \end{cases}\)          
dp_force \(\mathbf{F^{d}}\) Dashpot force (contact plane coord. system) VEC \(\mathbb{R}^3\) \(\mathbf{0}\) NO NO
    \(\left( -F_n^d,F_{ss}^d,F_{st}^d \right) \quad \left(\mbox{2D model: } F_{ss}^d \equiv 0 \right)\)          

Surface Property Inheritance

The shear modulus \(G\), Poisson’s ratio \(\nu\), and friction coefficient \(\mu\) may be inherited from the contacting pieces. Remember that for a property to be inherited, the inheritance flag for that property must be set to true (default), and both contacting pieces must hold a property with the exact same name.

The shear modulus and Poisson’s ratio are inherited as:

(18)\[\begin{split}\nu &= \frac{4 G^* - E^*}{2 G^* - E^*} \\ G &= 2 G^* (2 - \nu)\end{split}\]

with:

(19)\[\begin{split}E^* &= \left( \frac{1-\nu^{(1)}}{2G^{(1)}} + \frac{1-\nu^{(2)}}{2G^{(2)}} \right)^{-1} \\ G^* &= \left( \frac{2-\nu^{(1)}}{G^{(1)}} + \frac{2-\nu^{(2)}}{G^{(2)}} \right)^{-1}\end{split}\]

where (1) and (2) denote the properties of piece 1 and 2, respectively.

The friction coefficient is inherited using the minimum friction coefficient:

(20)\[\mu = \min(\mu^{(1)},\mu^{(2)})\]

Methods

No methods are defined by the Hertz model.

Callback Events


Table 3: Hertz Model Callback Events
Event Array Slot Value Type Range Description
contact_activated contact becomes active
  1 C_PNT N/A contact pointer
slip_change Slip state has changed
  1 C_PNT N/A contact pointer
  2 INT {0;1} slip change mode
        \(\;\;\;\;\;\;\begin{cases} \mbox{0: slip has initiated} \\ \mbox{1: slip has ended} \end{cases}\)

Usage and Verification Examples

The verification problem “Hertz Contact Model: Complex Loading Paths” exercises the Hertz contact model.

Revision 28 included a modification to the equations used to derive \(G\) and \(\nu\) from the piece properties. Users are encouraged to use the updated logic. However, the example below demonstrates how the default inheritance mechanism can be overridden to ensure reproducibility with previous models if needed.

model new
model domain extent -1 1
contact cmat default model hertz 
ball create id 1 rad 0.1 position-x -0.1
ball create id 2 rad 0.1 position-x  0.1
ball attribute density 2000.0
ball property 'hz_shear' 1.0e9 'hz_poiss' 0.3 'fric' 0.0 range id 1
ball property 'hz_shear' 2.0e9 'hz_poiss' 0.2 'fric' 0.5 range id 2
model clean
contact list property
model save 'default'
model new
model domain extent -1 1
contact cmat default model hertz inheritance hz_shear off hz_poiss off   
ball create id 1 rad 0.1 position-x -0.1
ball create id 2 rad 0.1 position-x  0.1
ball attribute density 2000.0
ball property 'hz_shear' 1.0e9 'hz_poiss' 0.3 'fric' 0.0 range id 1
ball property 'hz_shear' 2.0e9 'hz_poiss' 0.2 'fric' 0.5 range id 2

fish define catch_contact(c)
  local b1 = contact.end1(c)
  local b2 = contact.end2(c)
  contact.prop(c,'hz_shear') = 0.5*(ball.prop(b1,'hz_shear') + ...
                                   ball.prop(b2,'hz_shear'))
  contact.prop(c,'hz_poiss') = 0.5*(ball.prop(b1,'hz_poiss') + ...
                                   ball.prop(b2,'hz_poiss')) 
end
fish callback add catch_contact event contact_create
model clean
contact list property
model save 'overriden'
program return

In this example, the inheritance mechanism is deactivated for \(G\) and \(\nu\) in the contact model assignement table, and a FISH function is registered as a callback to the contact_create event. As a result, this function is executed any time a new contact is created in the model, and sets the contact model properties (see Contact Callback Events).

Model Summary

An alphabetical list of the hertz model properties is given here.

References

[Agnolin2007a]Agnolin, I., and J. N. Roux. “Internal states of model isotropic granular packings. I. Assembling process, geometry, and contact networks,” Phys. Rev., 76, 061302 (2007).
[Cundall1988a]Cundall, P. A. “Computer Simulations of Dense Sphere Assemblies,” Micromechanics of Granular Materials, pp. 113-123. M. Satake and J. T. Jenkins, eds. Amsterdam: Elsevier Science Publishers B.V. (1988).
[Ellata1996]Elata, D., and J. G. Berryman. “Contact force-displacement laws and the mechanical behavior of random packs of identical spheres,” Mechanics of Materials, 24, 229 (1996).
[Hunt1975]Hunt, K. H., and F. R. Crossley, “Coefficient of restitution interpreted as damping in vibroimpact,” J. Appl. Mech., 7, pp. 440-445 (1975).
[Mindlin1953](1, 2, 3) Mindlin, R. D., and H. Deresiewicz. “Elastic Spheres in Contact under Varying Oblique Forces,” J. Appl. Mech., 20, 327-344 (1953).
[Ramírez1999]Ramírez, R., Pöschel, T., Brilliantov, N. V., and T. Schwager, “Coefficient of restitution of colliding viscoelastic spheres,” Physical review E, vol. 60, no. 4, p. 4465, (1999).
[Tsuji1992]Tsuji, Y., Tanaka, T., and T. Ishida, “Lagrangian numerical simulation of plug flow of cohesionless particles in a horizontal pipe,” Powder technology, vol. 71, no. 3, pp. 239–250 (1992).