Hysteretic Model


The hysteretic contact model in PFC consists of a combination of the elastic portion of the Hertz model as described in the “Hertz Model” section, combined with an alternate dashpot group consisting of a nonlinear visco-elastic element in the normal direction. This model is referred to in commands and FISH as hysteretic.


Revision 28 included the following additional modification:

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 hysteretic model. Please refer to the documentation of the “Hertz Model” for further information.

Behavior Summary

The formulation is based on the literature review by [Machado2012]. Three different dashpot modes, which correspond to three different expressions for the dashpot force calculation, are provided.

Activity-Deletion Criteria

The hysteretic model defines an activity distance of 0.0. Contact activity is updated as:

(1)\[\begin{split}\begin{cases} g_c \leq 0.0 & \Rightarrow \mbox{contact active} \\ g_c > 0.0 & \Rightarrow \mbox{contact inactive} \end{cases}\end{split}\]

As discussed in the “Contact Resolution” section, force-displacement calculations are skipped for inactive contacts.

Force-Displacement Law

Similar to the Hertz model, the hysteretic model transmits a force and no moment. The force embeds a non-linear elastic component \(\mathbf{F^{h}}\) and a viscous component \(\mathbf{F^{d}}\)

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

The nonlinear Hertz force can be decomposed into a normal and a shear contribution:

(3)\[\mathbf{F^{h}} = -F_{n}^{h} {\kern 1pt} \hat{\mathbf{n}}_\mathbf{c} +\mathbf{F_{s}^{h}}\]

Its formulation is equivalent to that of the non-linear Hertz force described in the Hertz contact model.

The dashpot force has a normal component only:

(4)\[\mathbf{F^{d}} =-F_{n}^{d} {\kern 1pt} \hat{\mathbf{n}}_\mathbf{c}\]

where the formulation of \(F_{n}^{d}\) can be expressed in terms of a hysteresis damping factor \(\chi\) as:

(5)\[F_{n}^{d} = \chi \delta_{n}^{\alpha} \dot{\delta}_n\]

where \(\delta_{n} = -g_c\) represents physical overlap, \(\dot{\delta}_n = \mathbf{\dot{\delta}} {\kern 1pt} \hat{\mathbf{n}}_\mathbf{c}\) is the normal component of the relative translational velocity, and the hysteresis damping factor \(\chi\) depends on the dashpot mode \(M_d\), the target restitution coefficient \(e_n\), the Hertz normal force coefficient \(h_n\) (see here), and the normal impact velocity \(\dot{\delta}^{(-)}_n\) with:

(6)\[\begin{split}\chi = \begin{cases} \dfrac{1-e_n^2}{e_n} \frac{h_n}{\dot{\delta}^{(-)}_n} & \mbox{, } M_d = 0 \\ \dfrac{3\left(1-e_n^2\right)}{4} \frac{h_n}{\dot{\delta}^{(-)}_n} & \mbox{, } M_d = 1 \\ \dfrac{3\left(1-e_n^2\right)e^{2\left(1-e_n\right)}}{4} \frac{h_n}{\dot{\delta}^{(-)}_n} & \mbox{, } M_d = 2 \end{cases}\end{split}\]

Note that for the case \(M_d=1\), the target restitution coefficient used in the equation above is the maximal value between \(e_n\) and a minimal value of \(e^*_n\).

Energy Partitions

The hysteretic 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:

    (7)\[E_{k} = \frac{\alpha}{(\alpha+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:

    (8)\[\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:

    (9)\[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: Hysteretic 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


The properties defined by the hysteretic 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: Hysteretic model properties
Keyword Symbol Description Type Range Default Modifiable Inheritable
hysteretic Model name
Hertz Group:
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\) 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_mode \(M_d\) Dashpot mode [-] INT {0;1;2} 0 YES NO
    \(\;\;\;\;\;\;\begin{cases} \mbox{0: default mode} \\ \mbox{1: mode 1} \\ \mbox{2: mode 2} \end{cases}\)          
dp_en \(e_n\) Target restitution coef. [-] FLT \([0.0,1.0]\) 1.0 YES NO
dp_enmin \(e^*_n\) Minimal restitution coef. [-] FLT \([0.0,1.0]\) 0.0 YES NO
dp_force \(\mathbf{F^{d}}\) Dashpot force (contact plane coord. system) VEC \(\mathbb{R}^3\) \(\mathbf{0}\) NO NO
    \(\left( -F_n^d,0,0 \right)\)          

Property Inheritance

The shear modulus \(G\), Poisson’s ratio \(\nu\), and friction coefficient \(\mu\) may be inherited from the contacting pieces. The inheritance logic is similar to that of the hertz contact model described here.


No methods are defined by the hysteretic model.

Callback events

Table 3: Hysteretic Contact 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}\)

Model Summary

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

This is a hidden link to: [Gonthier2004]

And this is a hidden link to: [Lankarani1990]


[Gonthier2004] Y. Gonthier, J. McPhee, C. Lange, J-C. Piedboeuf, “A regularized contact model with asymmetric damping and dwell-time dependent friction”, Multibody System Dynamics, vol. 11, pp. 209-233 (2004).
[Lankarani1990]H.M. Lankarani, P.E. Nikravesh, “A contact force model with hysteresis damping for impact analysis of multibody systems”, Journal of Mechanical Design, vol. 112, pp. 369-376 (1990).
[Machado2012] M. Machado, P. Moreira, P. Flores, H.M. Lankarani, “Compliant contact force models in multibody dynamics: Evolution of the Hertz contact theory”, Mechanism and Machine Theory, vol. 53, pp. 99-121 (July 2012).