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 ballball and ballfacet 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 finitesized 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 forcedisplacement law is skipped for inactive contacts.
Note
The Hertz model implementation has changed from PFC 4.0 as follows:
BallFacet Contact Treatment: In PFC 4.0, the wall in a ballwall contact was conceptually replaced by an identical ball positioned symmetrically about the wall location so that the value of the overlap used in the forcedisplacement law was set equal to twice the physical ballwall overlap. In PFC 5.0 and PFC 6.0, the force is resolved for a ballfacet contact using the facet curvature of 0 (or, alternatively, the facet radius of \(\infty\)). This means that the normal force at a ballfacet contact is half the value calculated in PFC 4.0.
2D Formulation: The theoretical background of the Hertz model in PFC is only applicable to the case of 3D ballball and ballfacet contacts. Users may however use the Hertz contact model in PFC2D and select the value of the powerlaw exponent \(\alpha_h\) (see below).
The powerlaw 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.
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:
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.
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 nonzero value. Following the dimensional analysis of [Ramírez1999], \(\alpha_d = (\alpha_h1) / 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].
A new property, \(g_r\), has been added to allow the user to specify a nonzero reference gap.
ActivityDeletion 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,” forcedisplacement calculations are skipped for inactive contacts.
ForceDisplacement Law
The forcedisplacement law for the Hertz model updates the contact force and moment:
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:
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:
The forcedisplacement law for the Hertz model consists of the following steps:
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)\]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 Equation (12) 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_h1)/\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.
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 \; (notension\; 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_h1}}} \\ {\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 ballball}} \\ {\qquad m^{\left(1\right)} {\rm ,}\qquad {\rm ballfacet}} \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 nondimensionnal exponent, and \(\dot{\delta }_{n}\) is the relative normal translational velocity of this equation of the “Contact Resolution” section. For the dashpot normalbehavior mode of notension, \(F_{n}^{d}\) is capped to ensure that the total normal force (\(F_{n}^{l} +F_{n}^{d}\)) does not become tensile.
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 \; (slipcut)}} \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 ballball}} \\ {\qquad m^{\left(1\right)} {\rm ,}\qquad {\rm ballfacet}} \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 shearbehavior mode of slipcut, \(\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.
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}\]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 sheardisplacement 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.
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.
Keyword 
Symbol 
Description 
Range 
Accumulated 

Hertz Group: 


\(E_{k}\) 
strain energy 
\([0.0,+\infty)\) 
NO 

\(E_{\mu}\) 
total energy dissipated by slip 
\([0.0,+\infty)\) 
YES 
Dashpot Group: 


\(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.
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\) 
Shearforce 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:
with:
where (1) and (2) denote the properties of piece 1 and 2, respectively.
The friction coefficient is inherited using the minimum friction coefficient:
Methods
No methods are defined by the 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 positionx 0.1
ball create id 2 rad 0.1 positionx 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 positionx 0.1
ball create id 2 rad 0.1 positionx 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
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).
Cundall, P. A. “Computer Simulations of Dense Sphere Assemblies,” Micromechanics of Granular Materials, pp. 113123. M. Satake and J. T. Jenkins, eds. Amsterdam: Elsevier Science Publishers B.V. (1988).
Elata, D., and J. G. Berryman. “Contact forcedisplacement laws and the mechanical behavior of random packs of identical spheres,” Mechanics of Materials, 24, 229 (1996).
Hunt, K. H., and F. R. Crossley, “Coefficient of restitution interpreted as damping in vibroimpact,” J. Appl. Mech., 7, pp. 440445 (1975).
Mindlin, R. D., and H. Deresiewicz. “Elastic Spheres in Contact under Varying Oblique Forces,” J. Appl. Mech., 20, 327344 (1953).
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).
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).
Was this helpful? ...  Itasca Software © 2023, Itasca  Updated: Feb 16, 2024 