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 [Morrisey2013].

The implementation of this model in PFC is based on work from [Morrisey2013]. 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 ([Morrisey2013]):

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

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 this equation 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,MrμrˉR(FEEPAnFmin)μrˉR(FEEPAnFmin)(Mr/Mr),otherwise.

where the same arguments used for the shear force are used to increase the limiting moment by the effects of the minimum force Fmin. The normal force FEEPAn 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 EEPA 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((FEEPAn)o+FEEPAn)Δδn+12((FEEPAs)o+FEEPAs)ΔδδkswithΔδδks=ΔδδsΔδδμs=(FEEPAs(FEEPAs)okts)

    where (FEEPAn)o and (FEEPAs)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 (Δδδ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((FEEPAs)o+FEEPAs)Δδδμs
  3. Incrementally update the dashpot energy:

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

    where ˙δδ is the relative translational velocity of this equation in the i 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: EEPA Model Energies
Keyword Symbol Description Range Accumulated
Hertz Group:
energy‑strain‑eepa 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 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 gr Reference gap [length] FLT R 0.0 YES NO
eepa_shear G Shear modulus [stress] FLT (0.0,+) 0.0 YES YES
eepa_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 γ
Effective surface adhesion energy
[energy/area]
FLT [0.0,+) 0.0 YES NO
overlap_max δmax Maximum overlap [length] FLT [0.0,+) 0.0 NO N/A
force_min fmin Minimum force [force] FLT [0.0,+) 1.5 YES NO
pull_off F0 Pull-off force [force] FLT (,0.0] 0.0 NO NO
plas_ratio λp Plasticity ratio [-] FLT (0.0,1.0) 0.5 YES NO
lu_exp m Loading-unloading branch exponent [-] FLT [1.0,+) 1.5 YES NO
adh_exp χ Adhesive branch exponent [-] FLT [1.0,+) 1.5 YES NO
eepa_slip s Slip state [-] BOOL {false,true} false NO N/A
eepa_force FEEPA EEPA force (contact plane coord. system) VEC R3 0 YES NO
    (FEEPAn,FEEPAss,FEEPAst)(2D model: FEEPAss0)          
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 an identical 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(μ(1),μ(2))μr=min(μ(1)r,μ(2)r)

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

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.
[Marshall2009]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.
[Morrisey2013](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).