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], [Morrisey2013b], [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 [Marshall2009b]:

(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) ([Marshall2009b]):

(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:

(19)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:

(20)Mr:=MrktrΔθθb

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

The tangent rolling resistance stiffness ktr is defined as:

(21)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:

(22)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:

    (23)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:

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

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

    where ˙δδ is the relative translational velocity of Equation (10) of the “Contact Resolution” section.

  4. Update the rolling strain energy:

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

    (27)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:

(28)ν=4GE2GEG=2G(2ν)

with:

(29)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:

(30)μ=min

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

\;\;\;\;\;\;\begin{cases} \mbox{0: slip has initiated} \\ \mbox{1: slip has ended} \end{cases}

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.

[Marshall2009b] (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.

[Morrisey2013b]

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.