Edinburgh-Elasto-Plastic-Adhesive (EEPA) 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. Material subjected to relatively high loads (consolidation) may experience plastic deformation that results in increased contact area and increased cohesive behavior. History-dependent elasto-plastic-adhesive contact models can be used in DEM to simulate this behavior. One such model is the Edinburgh-Elasto-Plastic-Adhesion (EEPA) contact model developed by [Morrisey2013a].

The implementation of this model in PFC is based on work from [Morrisey2013a]. It was first proposed as a C++ Plugin for PFC 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 eepa.

Behavior Summary

The Edinburgh-Elasto-Plastic-Adhesion (EEPA) contact model is an extension of the linear hysteretic model by [WaltonBraun1986]. This model allows tensile forces to develop, as well as a non-linear force-displacement behavior in compression. The model also incorporates viscous damping and rolling resistance mechanisms similar to the Rolling Resistance Linear Model.

Activity-Deletion Criteria

Similar to the linear model, a contact with the EEPA model is active if and only if the surface gap is less than or equal to zero, and the force-displacement law is skipped for inactive contacts. The surface gap is shown in a figure in the linear model description. When the reference gap is zero, the notional surfaces coincide with the piece surfaces.

Force-Displacement Law

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

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

where FEEPA 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 FEEPA=FEEPAnˆnc+FEEPAs), the dashpot force, and the rolling resistance moment are updated as detailed below.

Normal Force FEEPAn

The model is an extension of the linear hysteretic model by [WaltonBraun1986], which allows tensile force development and a non-linear force-displacement behavior as shown in Figure 1 below.

The force-displacement law in the normal direction is defined by the pull-off force F0, the loading branch stiffness k1, the loading-unloading branch stiffness k2, the minimum force Fmin, the adhesion branch stiffness ka, and the plastic overlap δp.

../../../../../_images/cmeepa_fn.png

Figure 1: Force-displacement law in the normal direction.

The normal force is given by:

(2)FEEPAn={F0+k1δmnif k2(δmnδmp)k1δmnF0+k2(δmnδmp)if k1δmn>k2(δmnδmp)>kaδχnF0kaδχnδmnif kaδχnk2(δmnδmp)

The contact is activated at normal overlap δn=0, where the normal force jumps to a (negative) value F0. Virgin loading follows the k1-branch defined by the stiffness k1 and the overlap δn raised to the power m (the first condition in equation (2)). The maximum overlap is a history-dependent parameter and is updated and stored in the contact model.

Unloading follows the k2-branch defined by the stiffness k2 and the plastic overlap δp (the second condition in equation (2)). This branch uses the same exponent m as the k1-branch, and the stiffness k2 is defined in terms of the plasticity ratio λp, as

(3)k2=k11λp

When λp=0, then k2=k1 and there is no plastic behavior. When λp1, then k2+ and the behavior is perfectly plastic. The stiffness k1 is defined in terms of the effective Young’s modulus E, in a format similar to that of the classic Hertz model:

(4)k1=43EˉR

where

(5)ˉR=R1R2R1+R2

with R1 and R2 the radii of the two contacting pieces. If one piece is a wall facet, ˉR=Rp, where Rp is the particle radius.

The plastic overlap is defined as the overlap where the force on the k2-branch is equal to the pull-off force F0, and is given in terms of the maximum overlap and plasticity ratio,

(6)δp=λ1/mpδmax

Unloading along the k2-branch results in a tensile (adhesive) force development. At the plastic overlap δp the force is equal to the pull-off force F0. Further unloading follows the same branch until the magnitude of the tensile force is equal to the minimum force Fmin, which is given by:

(7)Fmin=F032πγa

where γ is the contact effective surface adhesion energy based on the Johnson-Kendall-Roberts (JKR) model for adhesive forces, which takes van-der-Waals effects into account([Johnson1971a]). The contact patch radius (assuming a circular contact area) a is given by Hertz’s theory, using the plastic overlap δp,

(8)a=2δpˉR

If the surface energy is set to zero, the adhesion branch is a horizontal line and Fmin=F0.

Once the force is equal to Fmin during unloading on the k2-branch, further unloading switches to the ka- or adhesion branch (third condition in equation (2)). The overlap exponent is χ and the adhesion stiffness is calculated as follows, to ensure that the unloading branch results in the force being equal to Fmin when δn=0:

(9)ka=FminF0δχmin

where the overlap corresponding to the minimum force is given by:

(10)δmin=(δmp+FminF0k2)1/m

Note that the adhesion stiffness ka is not constant, and is indirectly dependent on the load history (dependent on Fmin, which depends on the plastic overlap δp, which again depends on the maximum overlap δmax).

If the contact is unloading along the k2-branch, and then starts to reload, the same k2-branch is followed until the k1-branch is reached. Further virgin loading follows the k1-branch and will result in an increased and updated δmax (and δp). If the contact is unloading along the ka-branch, and then starts to reload, a k2-branch is followed as shown in Figure 1. In this case, the value of δp is updated to the point that this k2-branch is equal to the pull-off force F0.

On the k2-branch, the force has a minimum limit if unloading is allowed to continue until zero overlap δn=0, and, under some circumstances, it is possible the unloading branch would not reach the minimum force Fmin. The minimum force limit that can be achieved when unloading on the k2-branch is given by:

(11)Flimitmin=F0k2δmp

If the minimum force given by equation (7) is less than this limit (larger in magnitude) FminFlimitmin, the adhesion branch would never be reached and the solution for the minimum overlap δmin in equation (10) would have no real solution. This might occur when, for example, a relatively high surface energy is specified, and a specific contact has a relatively small maximum overlap. In order to avoid such a situation, the minimum force Fmin is modified and set to the average of the pull-off force F0 and limit force Flimitmin,

(12)Fmin={F032πγaif Fmin>Flimitmin12(F0+Flimitmin)otherwise

Shear Force FEEPAs

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

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

where (FEEPAs)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 a function of the normal overlap δn ([Morrisey2013a]):

(14)kts=ksf8GReδn

where G is the effective shear modulus and ksf is a scaling factor.

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 marks 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 (FHn+2Fpo), where FHn is the elastic force according to Hertz’s theory and Fpo is the pull-off (maximum tensile) force ([Marshall2009a]).

A similar approach is used in the EEPA model, and the updated force from equation (13) is limited by the following criteria:

(15)FEEPAs={FEEPAsif |FEEPAs|μ(FEEPAnFmin)μ(FEEPAnFmin)(FEEPAs|FEEPAs|)otherwise

where the normal force FEEPAn is given by equation (2), and the minimum force Fmin is given by equation (12). Note that Fmin would have a negative value, and will thus always increase the reference force (FEEPAnFmin) which can never be smaller than zero.

If the first condition in equation (15) 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={mk1δm1nif on the k1-branchmk2δm1nif on the k2-branchχkaδχ1nif on the ka-branch

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) in the i 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 (5)).

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

(22)Mr={Mr,

where the same arguments used for the shear force are used to increase the limiting moment by the effects of the minimum force F_{min}. The normal force F_n^{EEPA} is taken at the end of the current step, and \mu_r is the rolling coefficient of friction. If the first condition in (22) is met, then the rolling slip state is set to s_r=\mbox{false}, and when the second condition is met it is set to s_r=\mbox{true}.

Energy Partitions

The EEPA model provides five energy partitions:

  • strain energy, E_{k}, stored in the non-linear springs;

  • slip energy, E_{\mu}, defined as the total energy dissipated by frictional slip;

  • dashpot energy, E_{\beta}, defined as the total energy dissipated by the dashpots;

  • rolling strain energy, E_{k_r}, stored in the rolling linear spring; and

  • rolling slip energy, E_{\mu_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)\begin{split}\begin{array}{ll} E_k := E_k & + \frac{1}{2} \left(\left(F_{n}^{EEPA}\right)_o + F_{n}^{EEPA}\right) \Delta \delta_n \\ & + \frac{1}{2} \left(\left(\mathbf{F_{s}^{EEPA}}\right)_o + \mathbf{F_{s}^{EEPA}}\right) \cdot \Delta \pmb{δ}_\mathbf{s}^\mathbf{k} \\ {\rm with} & \Delta \pmb{δ} _\mathbf{s}^{k } =\Delta \pmb{δ} _\mathbf{s} -\Delta \pmb{δ} _\mathbf{s}^{\mu} = \left(\frac{\mathbf{F_{s}^{EEPA}} -\left(\mathbf{F_{s}^{EEPA}} \right)_{o} }{k_{s}^t } \right) \end{array}\end{split}

    where \left(F_{n}^{EEPA}\right)_o and \left(\mathbf{F_{s}^{EEPA}}\right)_o are, respectively, the EEPA normal and shear forces at the beginning of the step, and the relative shear-displacement increment has been decomposed into an elastic \left(\Delta \pmb{δ} _\mathbf{s}^{k} \right) and a slip \left(\Delta \pmb{δ} _\mathbf{s}^{\mu } \right) component, with k_s^t the tangent shear stiffness at the start of the current step (equation (14)).

  2. Incrementally update the slip energy:

    (24)E_{\mu } :=E_{\mu } - {\tfrac{1}{2}} \left(\left(\mathbf{F_{s}^{EEPA}} \right)_{o} +\mathbf{F_{s}^{EEPA}} \right)\cdot \Delta \pmb{δ} _\mathbf{s}^{\mu }
  3. Incrementally update the dashpot energy:

    (25)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 Equation (10) in the i Contact Resolution section.

  4. Update the rolling strain energy:

    (26)E_{k_r} = \frac{1}{2} \frac{\left\| \mathbf{M^{r}} \right\| ^{2} }{k_{r}}.
  5. Incrementally update the rolling slip energy:

    (27)\begin{split}\begin{array}{l} {E_{\mu_r} := E_{\mu_r} - {\tfrac{1}{2}} \biggl(\left(\mathbf{M^{r}} \right)_{o} +\mathbf{M^{r}} \biggr)\cdot \Delta \pmb{θ}_\mathbf{b}^{\mu_r} } \\ {{\rm with}\qquad \Delta \pmb{θ}_\mathbf{b}^{\mu_r} = \Delta \pmb{θ}_\mathbf{b} -\Delta \pmb{θ}_\mathbf{b}^{k} =\Delta \pmb{θ}_\mathbf{b} -\left(\frac{\mathbf{M^{r}} -\left(\mathbf{M^{r}} \right)_{o} }{k_{r}^t } \right)} \end{array}\end{split}

    where \left(\mathbf{M^{r}} \right)_{o} is the rolling resistance moment at the beginning of the timestep, and the relative bend-rotation increment has been decomposed into an elastic \Delta \pmb{θ}_\mathbf{b}^{k} and a slip \Delta \pmb{θ}_\mathbf{b}^{\mu_r} component, with k_{r}^t 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: EEPA Model Energies

Keyword

Symbol

Description

Range

Accumulated

Hertz Group:

energy‑strain‑eepa

E_{k}

strain energy

(-\infty,+\infty)

YES

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

Rolling-Resistance Group:

energy‑rrstrain

E_{k_r}

Rolling strain energy

[0.0,+\infty)

NO

energy‑rrslip

E_{\mu_r}

Total energy dissipated by rolling slip

(-\infty,0.0]

YES

Properties

The properties defined by the EEPA contact model are listed in the table below for a concise reference; see the “Contact Properties” section of the “Contact Model Framework” 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: EEPA Model Properties

Keyword

Symbol

Description

Type

Range

Default

Modifiable

Inheritable

eepa

Model name

EEPA Group:

rgap

g_r

Reference gap [length]

FLT

\mathbb{R}

0.0

YES

NO

eepa_shear

G

Shear modulus [stress]

FLT

(0.0,+\infty)

0.0

YES

YES

eepa_poiss

\nu

Poisson’s ratio [-]

FLT

[0.0,0.5]

0.0

YES

YES

ks_fac

k_{sf}

Shear stiffness scaling factor [-]

FLT

(0.0,+\infty)

1.0

YES

NO

fric

\mu

Sliding friction coefficient [-]

FLT

[0.0,+\infty)

0.0

YES

YES

surf_adh

\gamma^*

Effective surface adhesion energy

[energy/area]

FLT

[0.0,+\infty)

0.0

YES

NO

overlap_max

\delta_{max}

Maximum overlap [length]

FLT

[0.0,+\infty)

0.0

NO

N/A

force_min

f_{min}

Minimum force [force]

FLT

[0.0,+\infty)

1.5

YES

NO

pull_off

F^0

Pull-off force [force]

FLT

(-\infty,0.0]

0.0

NO

NO

plas_ratio

\lambda_p

Plasticity ratio [-]

FLT

(0.0,1.0)

0.5

YES

NO

lu_exp

m

Loading-unloading branch exponent [-]

FLT

[1.0,+\infty)

1.5

YES

NO

adh_exp

\chi

Adhesive branch exponent [-]

FLT

[1.0,+\infty)

1.5

YES

NO

eepa_slip

s

Slip state [-]

BOOL

{false,true}

false

NO

N/A

eepa_force

\mathbf{F^{EEPA}}

EEPA force (contact plane coord. system)

VEC

\mathbb{R}^3

\mathbf{0}

YES

NO

\left( -F_n^{EEPA},F_{ss}^{EEPA},F_{st}^{EEPA} \right) \quad \left(\mbox{2D model: } F_{ss}^{EEPA} \equiv 0 \right)

Rolling-Resistance Group:

rr_fric

\mu_r

Rolling friction coefficient [-]

FLT

[0.0,+\infty

0.0

YES

YES

rr_moment

M^r

Rolling Resistance moment (contact plane coord. system)

VEC

\mathbb{R}^3

\mathbf{0}

YES

NO

\left( 0,M_{bs}^r,M_{bt}^r \right) \quad \left(\mbox{2D model: } M_{bt}^r \equiv 0 \right)

rr_slip

s_r

Rolling slip state [-]

BOOL

{false,true}

false

NO

N/A

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_mode

M_d

Dashpot mode [-]

INT

{0;1}

0

YES

NO

\;\;\;\;\;\;\begin{cases} \mbox{0: no cutt-off} \\ \mbox{1: cut-off in shear when sliding} \end{cases}

dp_force

\mathbf{F^{d}}

Dashpot force (contact plane coord. system)

VEC

\mathbb{R}^3

\mathbf{0}

NO

N/A

\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, friction coefficient \mu and rolling friction coefficient \mu_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 an identical name.

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

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

with:

(29)\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 and rolling friction coefficients are inherited using the minimum of the values set for the pieces:

(30)\begin{split}\mu &= \min \left( \mu^{(1)} ,\mu^{(2)} \right) \\ \mu_r &= \min \left( \mu_r^{(1)},\mu_r^{(2)}\right)\end{split}

Methods

No methods are defined by the EEPA model.

Callback Events

Table 3: EEPA 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 EEPA model properties is given here.

References

[Johnson1971a] (1,2)

Johnson, K. L., Kendall, K., 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, 234(1558):301-313, Sept. 1971.

[Marshall2009a]

Marshall, J. S., “Discrete-Element Modeling of Particulate Aerosol Flows,” Journal of Computational Physics, 228(5):1541-1561, march 2009. doi:10.1016/j.jcp.2008.10.035.

[Morrisey2013a] (1,2,3)

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).

[Thornton1991]

Thornton, C., “Interparticle sliding in the presence of adhesion,” Journal of Physics D: Applied Physics, 24(11):1942-1946,1991. ISSN 13616463. doi:10.1088/0022-3727/24/11/007.

[ThorntonAndYin1991]

Thornton, C., and K. K. Yin, “Impact of elastic spheres with and without adhesion,” Powder Technology, 65(1-3):153-166,1991. ISSN 00325910. doi:10.1016/0032-5910(91)80178-L.

[WaltonBraun1986] (1,2)

Walton, R., and R. L. Braun. “Viscosity, Granular-Temperature, and Stress Calculations for Shearing Assemblies of Inelastic, Frictional Disks,” Journal of Rheology, 30(5), pp. 949-980 (1986).