Johnson-Kendall-Roberts (JKR) Contact Model

Introduction

Cohesion or adhesion can be modeled by introducing an attractive force component to the contact model. The terms cohesion and adhesion are usually defined for attracting forces between similar and dissimilar materials respectively. However, here the two terms are used interchangeably.

The Johnson-Kendall-Roberts (JKR) contact model is an extension of the well-known Hertz contact model proposed by [Johnson1971]. The model accounts for attraction forces due to van der Waals effects. The model is, however, also used to model material where the adhesion is caused by capillary or liquid-bridge forces ([Hærvig2017], [Carr2016],[Morrissey2013]_, [Xia2019]).

The implementation of this model in PFC was first proposed as a C++ Plugin for PFC version 6.0 by Prof. Corné Coetzee, whose contribution is gratefully acknowledged. The model for PFC 6.0 was made publicly available in the UDM Library on the Itasca website.

This model is incorporated as a built-in model in PFC 7.0. It is referred to in commands and FISH by the name jkr.

Behavior Summary

The JKR contact model is an extension of the Hertz-Mindlin model, which allows for tensile forces to develop due to surface adhesion. This model also incorporates viscous damping and rolling resistance mechanisms similar to those of the Rolling Resistance Linear Model.

Activity-Deletion Criteria

A contact with the JKR model becomes active when the surface gap becomes less or equal to the reference gap. Once activated, the contact remains active until the surface gap reaches a threshold positive tear-off distance. Note that a simplified version of the model may be used where the tear-off distance is set equal to the contact reference gap.

Force-Displacement Law

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

(1)Fc=FJKR+Fd,Mc=Mr

where FJKR is a non-linear force, Fd is the dashpot force, and Mr is the rolling resistance moment. The non-linear force (resolved into normal and shear components FJKR=FJKRnˆnc+FJKRs), the dashpot force, and the rolling resistance moment are updated as detailed below.

Normal Force FJKRn

Following [Johnson1971] and [Chokshi1993], an attractive force component can be introduced to the classic Hertz model, which acts over the circular contact patch area with radius a,

(2)Fadh=πa2γ

where γ is an effective surface energy per unit area. For two dissimilar surfaces the effective surface energy is given by

(3)γ=γ1+γ22γ12

where γ1 and γ2 are the surface energy of each surface, and γ12 is the interface energy. When the two surface materials are identical, γ1=γ2=γ, the interface energy is zero, γ12=0, and the effective surface energy reduces to

(4)γ=2γ

and the attractive force becomes:

(5)Fadh=2πa2γ

Due to the contact stress distribution the total contact normal force is given by:

(6)FJKRn=4Ea33ˉR16πγEa3

where the first term is identical to the elastic Hertz force (with ˉR the effective contact radius), and the second term is the adhesion component. Due to the attractive force, the contact area is larger than that predicted by the Hertz theory and the contact patch radius a is obtained by solving for the positive root of equation (6)

(7)a3=3ˉR4E[FJKRn+6πγˉR+FJKRn12πγˉR+(6πγˉR)2]

The Hertz relation between the contact patch radius a and the normal deformation (a=ˉRδn) no longer holds. According to the JKR theory, the normal deformation is related to the contact radius according to the following non-linear expression,

(8)δn=a2ˉR4πγaE

In PFC the contact overlap δn is readily available, and from this the contact patch radius a should be calculated. Solving equation (8) for the contact patch radius is however not trivial. [Parteli2014] proposed an analytical solution based on a fourth-order expansion of this equation and solving for the real root that is larger than the patch radius of the classic Hertz model, which yields

(9)a=12(w+w24(c2+s+λ)withc0=ˉR2δ2nc1=8(1ν2)πγˉR2Ec2=2ˉRδnP=c2212c0Q=c32108+c0c23c218U=(Q2+Q24+P327)1/3s={5c26+UP3Uif P05c26Q1/3if P=0w=c2+2sλ=c12w

This value of a is then used to compute FJKRn using equation (6).

The maximum tensile force is defined as the pull-off force Fpo ([Chokshi1993]),

(10)Fpo=3πγˉR

The pull-off force is independent on the elastic properties and only a function of particle size and surface adhesion energy. The contact force FJKRn can be normalized using Fpo and plotted as a function of the normalized overlap δn/δto, where the tear-off distance δto is the separation distance at which the contact breaks under tension:

(11)δto=12161/3a20ˉR

where a0 is the contact patch radius where FJKRn=0, i.e., at the point where the contact is in equilibrium when no other external forces are acting:

(12)a0=(9πγˉR2E)1/3

Plotting the normalized force FJKRn/Fpo an normalized patch radius a/a0 against the normalized overlap δn/δto results in Figure 1 below. When two particles approach each other, the contact is formed when the two particles physically touch at δn=0 (point 3 in Figure 1). At this point, the force suddenly jumps to a value FJKRn=8/9Fpo, with the contact patch radius a=(2/3)2/3a0. This attractive force will pull the particles closer until the equilibrium point is reached (assuming no other forces are acting) where FJKRn=0, a=a0, and δn=(4/3)2/3δto (point 4 in Figure 1). Further loading due to external forces will result in an increased positive overlap δn and an increase in the compressive force as defined by equation (6).

Upon unloading, the same path is followed up to zero overlap δn=0, after which the force will further decrease down to a minimum value FJKRn=Fpo due to the necking effect where a=(1/2)2/3a0 (point 2 in Figure 1). With further unloading, the force increases to a value FJKRn=5/9Fpo where the contact suddenly breaks (tear-off) and the force jumps to zero (point 1 in Figure 1). At this point the overlap is given by δn=δto and the patch radius by a=(1/6)2/3a0.

../../../../../_images/cmjkr_fn.png

Figure 1: Normalized force and normalized contact patch radius of the JKR model as a function of the normalized overlap.

Note that the normal force is in tension when the overlap is positive δn>0 and between points 3 and 4 in Figure 1. This is different from the Hertz model, where the normal force is always compressive when the overlap is positive.

Also, in the Hertz model, the normal force becomes zero for δn0, while in the JKR model the contact remains active with a finite contact patch until the tear-off separation distance is reached (point 1 in Figure 1).

A simplified version of the model can be used where the normal force becomes zero for δn0 by setting Ma=0.

Shear Force FJKRs

In the shear direction, a trial shear force is first computed as:

(13)FJKRs=(FJKRs)o+ktsΔδδs

where (FJKRs)o is the shear force at the beginning of the timestep, Δδδs is the relative shear-displacement increment, and the tangent shear stiffness kts is given by [Marshall2009]:

(14)kts=8Ga

where G is the effective shear modulus. The tangent shear stiffness is a function of the patch radius a and not the normal overlap as in the Hertz model, because the patch radius can be present for negative overlap until the contact breaks.

To account for cohesive effects in the tangential direction, the Coulomb friction limit is usually modified. This is a relatively simple approach that is easy to implement and is based on the work by [Thornton1991] and [ThorntonAndYin1991], who showed that, with adhesion present, the contact area decreases with an increase in the tangential force. The tangential force reaches a critical value, which defines the transition from a “peeling” action to a “sliding” action. Assuming the JKR theory ([Johnson1971a]), it is shown that at this critical point the contact area corresponds to a Hertzian-like normal stress distribution under a load of (FJKRn+2Fpo) ([Marshall2009]):

(15)FJKRs={FJKRs,if |FJKRs|μ(FJKRn+2Fpo)μ(FJKRn+2Fpo)(FJKRs|FJKRs|),otherwise

If the first condition is met, the boolean slip state is set to s=false, and when the second condition is met it is set to s=true. Whenever the slip states changes, the slip_change callback event occurs.

Dashpot force

The dashpot force is resolved into normal and shear components:

(16)Fd=Fdnˆnc+Fds

In the normal direction:

(17)Fdn=2βnmcktn˙δn

where mc=m1m2m1+m2 is the effective contact mass with m1 and m2 the mass of the two contacting particles (for a contact with a wall facet, mc=mp where mp is the mass of the particle), ˙δn is the relative normal translational velocity, ktn is the tangent normal stiffness, and βn is the normal critical damping ratio. This ratio takes a value βn=0 when there is no damping, and βn=1 when the system is critically damped.

The tangent normal stiffness is given by:

(18)ktn=2Ea

The dashpot force in the shear direction is given by:

Fds={(2βsmckts)˙δδs,s= false or Md=0 (full shear)0,s= true and Md=1 (slip-cut)

where Md is the dashpot behavior mode, ˙δδs is the relative shear translational velocity, βs is the shear critical damping ratio, and kts is the tangential shear stiffness given by equation (14).

Rolling Resistance

The rolling resistance moment is first incremented as:

(19)Mr:=MrktrΔθθb

where ktr is the tangent rolling resistance stiffness, and Δθθb is the relative bend-rotation increment of this equation of the “Contact Resolution” section.

The tangent rolling resistance stiffness ktr is defined as:

(20)ktr=ktsˉR2

with kts the tangent shear stiffness given by equation (14) and ˉR the contact effective radius(equation :eq:’cmjkr_eqR’) CS: maybe this is cmeepa_eqR??? see line 262 (the cmeepa_Mr2 equation reference) for more badness).

The magnitude of the updated rolling resistance moment is then checked against a threshold limit:

(21)Mr={Mr,if MrμrˉR(FJKRn+2Fpo)μrˉR(FJKRn+2Fpo)(MrMr),otherwise.

where the same arguments as used for the shear force are used to increase the limiting moment by the effects of the pull-off force Fpo. The normal force FJKRn is taken at the end of the current step, and μr is the rolling coefficient of friction. If the first condition in (22) is met, then the rolling slip state is set to sr=false, and when the second condition is met it is set to sr=true

Energy Partitions

The JKR model provides five energy partitions:

  • strain energy, Ek, stored in the non-linear springs;
  • slip energy, Eμ, defined as the total energy dissipated by frictional slip;
  • dashpot energy, Eβ, defined as the total energy dissipated by the dashpots;
  • rolling strain energy, Ekr, stored in the rolling linear spring; and
  • rolling slip energy, Eμr, defined as the total energy dissipated by rolling slip.

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

  1. Incrementally update the strain energy:

    (22)Ek:=Ek+12((FJKRn)o+FJKRn)Δδn+12((FJKRs)o+FJKRs)ΔδδkswithΔδδks=ΔδδsΔδδμs=(FJKRs(FJKRs)okts)

    where (FJKRn)o and (FJKRs)o are, respectively, the JKR normal and shear forces at the beginning of the step, and the relative shear-displacement increment has been decomposed into an elastic (Δδδks) and a slip (Δδδμs) component, with kts the tangent shear stiffness at the start of the current step (equation (14)).

  2. Incrementally update the slip energy:

    (23)Eμ:=Eμ12((FJKRs)o+FJKRs)Δδδμs
  3. Incrementally update the dashpot energy:

    (24)Eβ:=EβFd(˙δδΔt)

    where ˙δδ is the relative translational velocity of this equation of the “Contact Resolution” section.

  4. Update the rolling strain energy:

    (25)Ekr=12Mr2kr.
  5. Incrementally update the rolling slip energy:

    (26)Eμr:=Eμr12((Mr)o+Mr)ΔθθμrbwithΔθθμrb=ΔθθbΔθθkb=Δθθb(Mr(Mr)oktr)

    where (Mr)o is the rolling resistance moment at the beginning of the timestep, and the relative bend-rotation increment has been decomposed into an elastic Δθθkb and a slip Δθθμrb component, with ktr the tangent rolling stiffness. The dot product of the slip component and the rolling resistance moment occurring during the timestep gives the increment of rolling slip energy.

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

Table 1: JKR Model Energies
Keyword Symbol Description Range Accumulated
JKR Group:
energy‑strain‑jkr Ek strain energy (,+) YES
energy‑slip Eμ total energy dissipated by slip [0.0,+) YES
Dashpot Group:
energy‑dashpot Eβ total energy dissipated by dashpots [0.0,+) YES
Rolling-Resistance Group:
energy‑rrstrain Ekr Rolling strain energy [0.0,+) NO
energy‑rrslip Eμr Total energy dissipated by rolling slip (,0.0] YES

Properties

The properties defined by the JKR contact model are listed in the table below for 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: JKR Model Properties
Keyword Symbol Description Type Range Default Modifiable Inheritable
jkr Model name
JKR Group:
jkr_shear G Shear modulus [stress] FLT (0.0,+) 0.0 YES YES
jkr_poiss ν Poisson’s ratio [-] FLT [0.0,0.5] 0.0 YES YES
ks_fac ksf Shear stiffness scaling factor [-] FLT (0.0,+) 1.0 YES NO
fric μ Sliding friction coefficient [-] FLT [0.0,+) 0.0 YES YES
surf_adh γ Surface adhesion energy [energy/area] FLT [0.0,+) 0.0 YES NO
a0 a0 Equilibrium contact patch radius [length] FLT [0.0,+) 0.0 NO N/A
pull_off_f Fpo Pull-off force [force] FLT (,0.0] 0.0 NO NO
tear_off_d δto Tear-off distance [-] FLT (0.0,1.0) 0.5 YES NO
active_mode Ma Active mode [-] INT {0;1} 0 YES NO
    {0: no negative overlap1: allows negative overlap          
jkr_slip s Slip state [-] BOOL {false,true} false NO N/A
jkr_force FJKR JKR force (contact plane coord. system) VEC R3 0 YES NO
    (FJKRn,FJKRss,FJKRst)(2D model: FJKRss0)          
Rolling-Resistance Group:
rr_fric μr Rolling friction coefficient [-] FLT [0.0,+ 0.0 YES YES
rr_moment Mr Rolling Resistance moment (contact plane coord. system) VEC R3 0 YES NO
    (0,Mrbs,Mrbt)(2D model: Mrbt0)          
rr_slip sr Rolling slip state [-] BOOL {false,true} false NO N/A
Dashpot Group:
dp_nratio βn Normal critical damping ratio [-] FLT [0.0,1.0] 0.0 YES NO
dp_sratio βs Shear critical damping ratio [-] FLT [0.0,1.0] 0.0 YES NO
dp_mode Md Dashpot mode [-] INT {0;1} 0 YES NO
    {0: no cutt-off1: cut-off in shear when sliding          
dp_force Fd Dashpot force (contact plane coord. system) VEC R3 0 NO N/A
    (Fdn,Fdss,Fdst)(2D model: Fdss0)          

Surface Property Inheritance

The shear modulus G, Poisson’s ratio ν, friction coefficient μ, and rolling friction coefficient μr 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:

ν=4GE2GEG=2G(2ν)

with:

E=(1ν(1)2G(1)+1ν(2)2G(2))1G=(2ν(1)G(1)+2ν(2)G(2))1

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

The friction and rolling friction coefficients are inherited using the minimum of the values set for the pieces:

μ=min(μ(1),μ(2))μr=min(μ(1)r,μ(2)r)

Methods

No methods are defined by the JKR model.

Callback Events

Table 3: JKR 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
        {0: slip has initiated1: slip has ended

Model Summary

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

References

[Carr2016]Carr, M. J., W. Chen, K. Williams, and A. Katterfeld. “Comparative investigation on modelling wet and sticky material behaviours with a simplified JKR cohesion model and liquid bridging cohesion model in DEM,” ICBMH 2016 - 12th International Conference on Bulk Materials Storage, Handling and Transportation, Proceedings, pages 40–49, 2016.
[Chokshi1993](1, 2) Chokshi, A.G., Tielens, G. M., and D. Hollenbach. “Dust Coagulation,” The Astrophysical Journal, 407:806-819, Apr. 1993.
[Hærvig2017]Hærvig, J. , U. Kleinhans, C.Wieland, H. Spliethoff, A. L. Jensen, K. Sørensen, and T. J. Condra. “On the adhesive JKR contact and rolling models for reduced particle stiffness discrete element simulations,” Powder Technology, 319:472–482, 2017. ISSN 1873328X. doi:10.1016/j.powtec.2017.07.006. URL http://dx.doi.org/10.1016/j.powtec.2017.07.006.
[Johnson1971](1, 2) Johnson, K. L., K. Kendall, and A. D. Roberts. “Surface Energy and the Contact of Elastic Solids,” Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, 324(1558):301–313, Sept. 1971.
[Marshall2009](1, 2) Marshall, J. S., “Discrete-Element Modeling of Particulate Aerosol Flows,” Journal of Computational Physics, 228(5):1541-1561, Mar. 2009, doi:10.1016/j.jcp.2008.10.035.
[Morrisey2013]Morrissey, J. P. “Discrete Element Modelling of Iron Ore Pellets to Include the Effects of Moisture and Fines,” PhD thesis, Edinburgh, Scotland: University of Edinburgh (2013).
[Parteli2014]Parteli, E. J., Schmidt, J., Blümel, C., Wirth, K. E., Peukert, W., and T. Pöshel. “Attractive particle interaction forces and packing density of fine glass powders,” Scientific Reports, 4:1-7,2014. ISSN 20454322. doi:10.1038/srep06227.
[Xia2019]Xia, R., B. Li, X.Wang, T. Li, and Z. Yang. “Measurement and calibration of the discrete element parameters of wet bulk coal,” Measurement: Journal of the International Measurement Confederation, 142:84–95, 2019. ISSN 02632241. doi: 10.1016/j.measurement.2019.04.069.