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:
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.
The normal force is given by:
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
When λp=0, then k2=k1 and there is no plastic behavior. When λp→1, 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:
where
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,
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:
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,
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:
where the overlap corresponding to the minimum force is given by:
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:
If the minimum force given by equation (7) is less than this limit (larger in magnitude) Fmin≤Flimitmin, 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,
Shear Force FEEPAs
In the shear direction, a trial shear force is first computed as:
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]):
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:
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 (FEEPAn−Fmin) 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:
In the normal direction:
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:
The dashpot force in the shear direction is given by:
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:
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:
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:
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:
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)).
Incrementally update the slip energy:
(24)Eμ:=Eμ−12((FEEPAs)o+FEEPAs)⋅ΔδδμsIncrementally update the dashpot energy:
(25)Eβ:=Eβ−Fd⋅(˙δδΔt)where ˙δδ is the relative translational velocity of this equation in the i Contact Resolution section.
Update the rolling strain energy:
(26)Ekr=12‖Mr‖2kr.Incrementally update the rolling slip energy:
(27)Eμr:=Eμr−12((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.
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.
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 | γ∗ |
|
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: FEEPAss≡0) | |||||||
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: Mrbt≡0) | |||||||
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: Fdss≡0) |
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:
with:
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:
Methods
No methods are defined by the 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). |
Was this helpful? ... | FLAC3D © 2019, Itasca | Updated: Feb 25, 2024 |