JKR Adhesion and Capillary Bridging Interaction Model for DEM

DEM Adhesion Implementation

This tutorial explains a Discrete Element Method (DEM) interaction model designed for simulating cohesive granular materials. It combines the Johnson-Kendall-Roberts (JKR) theory for adhesion, a Coulomb-limited friction model, and includes provisions for capillary forces (currently commented out in the provided code snippet).

Note

Assumptions & Scope:

• The model simulates interactions between pairs of spherical particles.

• It calculates normal forces based on JKR theory (elasticity + adhesion) and damping.

• It calculates tangential forces based on a friction model (detailed in Section 2e) limited by Coulomb’s law.

• Damping (normal and tangential) is included, based on the coefficient of restitution.

• A capillary force model avilable (but disabled). It is intended to handle attraction due to liquid bridges when particles are close but not in contact.

• An optional, simple electrostatic force (Coulomb’s Law) is included for non-contacting particles.

1. Theoretical Background

a) Hertz Contact Theory (Baseline)

Hertz theory describes the elastic contact between two non-adhesive spheres. It forms the basis for many DEM contact models.

• Contact Area: The radius of the circular contact area, \(a\), depends on the overlap \(\delta\) (how much the particles deform into each other) and the effective radius \(R^*\).

\[a = \sqrt{R^* \delta}\]

• Normal Force: The repulsive normal force \(F_{\text{Hertz}}\) is proportional to the overlap raised to the power of 3/2.

\[F_{\text{Hertz}} = \frac{4}{3} E^* \sqrt{R^*} \delta^{3/2}\]

where \(E^*\) is the effective Young’s Modulus.

Hertz theory assumes purely repulsive elastic forces and no adhesion.

b) Johnson-Kendall-Roberts (JKR) Theory

JKR theory extends Hertz theory by incorporating the effect of surface energy (\(\gamma\), denoted SurfaceEnergy in the code) which leads to adhesion.

• Addition from Hertz Contact Theory: Surface attraction forces pull the particles together, increasing the contact area compared to Hertz theory for the same load, and allowing the contact to sustain a tensile (negative) force before separation.

• Contact Radius (\(a\)): The relationship between overlap (\(\delta\)), contact radius (\(a\)), effective radius (\(R^*\)), effective modulus (\(E^*\)), and surface energy (\(\gamma\)) is implicit:

  \[\delta = \frac{a^2}{R^*} - \sqrt{\frac{8 \pi a \gamma}{E^*}}\]

  This equation must typically be solved numerically to find \(a\) for a given \(\delta\).

• Normal Force (\(F_{\text{JKR}}\)): The normal force includes both the elastic repulsion (Hertz-like term) and the adhesion term:

  \[F_{\text{JKR}} = \frac{4}{3} E^* \frac{a^3}{R^*} - \sqrt{6 \pi \gamma E^* a^3}\]

  or alternatively, expressed using the Hertzian force \(F_{\text{Hertz}}(a)\) based on contact radius \(a\):

  \[F_{\text{JKR}} = F_{\text{Hertz}}(a) - \sqrt{6 \pi \gamma E^* a^3}\]

• Pull-off Force: The maximum tensile force required to separate the particles is independent of the elastic properties and given by:

  \[F_{\text{pull-off}} = -\frac{3}{2} \pi \gamma R^*\]

c) Capillary Bridging Theory

When a small amount of liquid is present (e.g., due to humidity), liquid bridges can form between nearby particles, creating an attractive capillary force. This force is significant in cohesive granular materials.

• Origin: This force arises from the surface tension (\(\sigma\), SurfTension in the code) of the liquid and the negative Laplace pressure inside the curved meniscus of the liquid bridge.

• Force Magnitude: The magnitude depends complexly on the separation distance (\(s\)), the effective radius \(R^*\), surface tension \(\sigma\), the liquid-solid contact angle \(\theta\) (ContactAngle), and the volume of the liquid bridge (related to relative humidity or liquid content).

• Theoretical Models: Detailed analytical or semi-analytical expressions for the capillary force, often derived using the Laplace-Young equation or energy minimization approaches, can be found in literature (e.g., Willett et al., 2000). These models capture the force variation with separation distance. Simplified approximations are often used in DEM codes for computational efficiency. The model from Willett et al. is implementated in the code.

• Rupture Distance: The liquid bridge can only exist up to a certain critical separation distance (s_Critical in the code), beyond which it ruptures, and the capillary force drops abruptly to zero. This rupture distance also depends on the factors mentioned above.

2. Implementation Details

a) Parameters and Initial Calculations

• Material Properties:

  • youngsM: Young’s Modulus (\(E\)) [Pa]

  • pr: Poisson’s Ratio (\(\nu\)) [-]

  • SurfaceEnergy: Surface energy (\(\gamma\)) [N/m or J/m²]

  • rho_p: Particle density (\(\rho_p\)) [kg/m³]

  • SurfTension: Liquid surface tension (\(\sigma\)) [N/m]

  • ContactAngle: Liquid-solid contact angle (\(\theta\)) [radians]

• Contact Properties:

  • COR: Coefficient of restitution [-] (Used for damping)

  • mu: Coefficient of friction (\(\mu\)) [-]

• Derived Parameters:

  • E_star: Effective Young’s Modulus:

    \[E^* = \frac{E}{2 (1 - \nu^2)}\]

  • Reff: Effective Radius:

    \[R^* = \frac{R_1 R_2}{R_1 + R_2} \quad (\text{where } R = d/2)\]

  • Meff: Effective Mass:

    \[M^* = \frac{m_1 m_2}{m_1 + m_2}\]

• Geometric/Kinematic:

  • r[3]: Vector from center of particle 1 to particle 2.

  • relvel[3]: Relative velocity vector (\(\mathbf{v}_{\text{rel}} = \mathbf{v}_1 - \mathbf{v}_2\)).

  • Distance: Magnitude of r.

  • UnitVect[3]: Unit vector (\(\mathbf{n}\)) along r (pointing from 2 to 1 based on the code’s definition of r).

  • touch_distance: Sum of radii (\(R_1 + R_2\)).

  • overlap: \(\delta = \text{touch\_distance} - \text{Distance}\). Positive for contact, negative for separation \(s = -\delta\).

b) Contact Detection

Contact is detected if the overlap is positive: if (overlap > 0.f)

c) Contact Force Calculation (overlap > 0.f)

1. Calculate JKR Contact Radius (\(a\)):

   • The code uses a numerical root-finding method (bisection) to solve the implicit JKR equation \(\text{residual}(a, ...) = 0\). The residual function represents the standard JKR equation rearranged:

     \[\text{residual}(a) = \delta - \frac{a^2}{R^*} + \sqrt{\frac{8 \pi a \gamma}{E^*}} \quad (= 0 \text{ at solution})\]

   • The bisection function iteratively finds an \(a\) that makes residual close to zero (within tolerance). * (Mapping to code variables: \(\delta\) = overlap, \(R^*\) = Reff, \(\gamma\) = SurfaceEnergy, \(E^*\) = E_star)*.

   • Error handling is included.

   • If bisection returns NaN, it defaults to the Hertzian contact radius: \(a = \sqrt{R^* \delta}\).

   • An alternative analytical solution (a_return based on Parteli (2014)) is provided but is less stable and likely requires smaller time steps for stability.

2. Calculate Normal Stiffness (\(k_n\)):

   • A stiffness \(k_n\) is calculated based on the Hertzian-like relation using the JKR contact radius \(a\):

     \[k_n = \frac{4}{3} E^* a\]

     (Note: This is used in normal and tangential damping calculations).

3. Calculate Normal Damping Coefficient (\(C_n\)):

   • Damping opposes the relative normal velocity, related to COR.

   • Damping Ratio (\(\beta\)):

     \[\beta = \frac{-\ln(COR)}{\sqrt{\pi^2 + (\ln(COR))^2}}\]

   • Damping Coefficient (DampCoeff):

     \[C_n = 2 \beta \sqrt{M^* k_n}\]

4. Calculate Normal Force (Norm_Force):

   • Combines the JKR force (elastic repulsion + adhesion) and the normal damping force:

     \[F_n = \underbrace{\left( \frac{4}{3} E^* \frac{a^3}{R^*} - \sqrt{6 \pi \gamma E^* a^3} \right)}_{\text{JKR Force}} + \underbrace{\left( -C_n (\mathbf{v}_{\text{rel}} \cdot \mathbf{n}) \right)}_{\text{Normal Damping Force}}\]

d) Non-Contact Force Calculation (overlap <= 0.f)

1. Electrostatic Force (Attractive_Force): * Enabled: A flag ChargeEnabled must be true.

   • Calculates force based on Coulomb’s law:

     \[F_{\text{elec}} = -k \frac{q_1 q_2}{\text{Distance}^2}\]

e) Capillary Force Implementation (Willett Model)

This section details the implementation of the attractive capillary force based on the approach described by Willett et al. (2000), using a parameterized dimensionless force function \(F^*\).

Note

This implementation requires the dimensionless liquid bridge volume \(V^* = V / (R^*)^3\) (represented by the variable vStar in the code) as an input or a pre-calculated parameter, where \(V\) is the actual liquid bridge volume.

1. Conditions for Activation:

   • The capillary force calculation is triggered by an if statement checking several conditions:

     • Enabled: A flag CapillaryForceEnabled must be true.

     • Separation: Particles must be separated (\(s = -\text{overlap} > 0\)), exceeding a small numerical tolerance (e.g., 1e-15 in the code) to avoid issues at exact contact.

     • Proximity: The separation \(s\) must be less than a critical rupture distance \(s_{\text{critical}}\) (s_critical).

     • Liquid Volume: A non-negligible liquid bridge volume must be present (\(V^*\) or vStar must be greater than a small tolerance, e.g., 1e-12).

2. Calculate Dimensionless Parameters:

   • If the activation conditions are met, dimensionless variables central to the Willett et al. model are computed:

     • Dimensionless separation (\(s^*\)): The ratio of separation to effective radius.

       \[s^* = s / R^*\]

     • Scaled dimensionless separation (\(s^+\)): Dimensionless separation scaled by the square root of dimensionless volume.

       \[s^+ = s^* / \sqrt{V^*}\]

3. Core Force Calculation (using ``FStar`` Function):

   • The capillary force magnitude (\(F_{\text{capillary}}\)) is calculated by multiplying a standard scaling factor by a dimensionless force function, \(F^*\). This function encapsulates the complex physics dependent on geometry and material properties.

     \[F_{\text{capillary}} = 2 \pi R^* \sigma F^*(\theta, V^*, s^+)\]

   • The code relies on an external function FStar which takes the contact angle ContactAngle (\(\theta\)), dimensionless volume vStar (\(V^*\)), and scaled separation s_plus (\(s^+\)) as input and returns the dimensionless force \(F^*\).

   • Note: The implementation of the FStar function itself is crucial and must accurately represent the findings of Willett et al. (2000), typically through curve-fitting formulae (like the polynomials they provide for specific contact angles), look-up tables, or other numerical approximations based on their theoretical results.

4. Apply Stability Limit (Heuristic):

   • To mitigate potential numerical instability from large attractive forces over small separations, the calculated capillary force is capped by a heuristic limit, CapLimit.

   • This limit is estimated based on the force required to move the reduced mass (\(M^*\) or Meff) across a fraction of the gap distance (\(s = -\text{overlap}\)) within a single time step (\(dt\)). The formula used in the code is:

     \[F_{\text{limit, stab}} = \frac{separation M^*}{2 (dt)^2}\]

   • The capillary force is then limited:

     \[F_{\text{capillary}} = \min(F_{\text{capillary}}, F_{\text{limit, stab}})\]

5. Force Application:

   • The final, potentially capped, capillary force magnitude \(F_{\text{capillary}}\) is added to the total interaction force vector, acting as an attraction along the line connecting the centers (i.e., in the direction \(-\mathbf{n}\) or -UnitVect).

     • fx -= F_capillary * UnitVect[0]; etc. applied after capping

f) Tangential Force (Friction)

Friction acts to oppose the relative tangential motion between contacting particles. This implementation uses a non-history model, calculating the kinetic friction force directly based on the Coulomb limit and the instantaneous relative tangential velocity direction. It does not track elastic tangential displacement (stick-slip deformation).

Note

This model assumes that any non-zero relative tangential velocity results in fully developed kinetic friction opposing the motion. Static friction (stiction requiring a force threshold to initiate sliding) is not explicitly modeled by tracking elastic deformation.

1. Calculate Relative Tangential Velocity and Direction:

   • The code calculates the relative tangential velocity \(\mathbf{v}_{t, rel}\) (stored in relvel_t).

   • If \(|\mathbf{v}_{t, rel}|\) is greater than a small tolerance (to avoid division by zero), the unit vector \(\mathbf{t}\) in the direction of relative tangential motion is determined (stored in tangentialUnitVect):

     \[\mathbf{t} = \frac{\mathbf{v}_{t, rel}}{|\mathbf{v}_{t, rel}|}\]

2. Calculate Coulomb Friction Limit:

   • The magnitude of the kinetic friction force is determined by the Coulomb friction law, based on the coefficient of friction \(\mu\) (mu in the code) and the normal force.

   • The normal force used for the limit typically only considers the compressive component, ensuring adhesion doesn’t increase the friction limit (though different models exist). We use \(\max(0, F_n)\) where \(F_n\) is the total normal force calculated previously (Norm_Force).

     \[F_{limit} = \mu \max(0, F_n)\]

3. Determine Tangential Force Vector:

   • In this non-history model, the magnitude of the kinetic friction force \(| \mathbf{F}_t |\) is set directly to the Coulomb limit whenever relative tangential motion occurs (i.e., when \(|\mathbf{v}_{t, rel}|\) is non-zero).

     \[| \mathbf{F}_t | = F_{limit} = \mu \max(0, F_n) \quad (\text{if } |\mathbf{v}_{t, rel}| > 0)\]

   • If there is no relative tangential velocity (\(|\mathbf{v}_{t, rel}| = 0\)), the tangential force is zero.

   • The friction force vector \(\mathbf{F}_t\) opposes the direction of relative tangential velocity \(\mathbf{t}\).

     \[\mathbf{F}_t = - | \mathbf{F}_t | \mathbf{t}\]

   • The calculated magnitude \(| \mathbf{F}_t |\) is presumably stored in the Tangential_Force variable in the code before being used in the final force assembly (Section 2f).

f) Assemble Total Force Vector

The total interaction force \(\mathbf{F}\) applied to the particles is the vector sum of the normal force (\(F_n \mathbf{n}\)) and the tangential friction force (\(\mathbf{F}_t = -|F_t| \mathbf{t}\)). The resulting force components are often stored in variables like fx, fy, and fz.

• During Contact (overlap > 0):

  The force components are calculated by projecting the normal and tangential forces onto the Cartesian axes:

  \[\mathbf{F} = F_n \mathbf{n} + \mathbf{F}_t = F_n \mathbf{n} - |F_t| \mathbf{t}\]

  The components (fx, fy, fz) are: \(f_x = F_n n_x - |F_t| t_x\) \(f_y = F_n n_y - |F_t| t_y\) \(f_z = F_n n_z - |F_t| t_z\)

  Where \(F_n\) is Norm_Force, \(|F_t|\) is the calculated magnitude of the tangential force (presumably stored in Tangential_Force in the updated code), \(\mathbf{n}\) is UnitVect, and \(\mathbf{t}\) is tangentialUnitVect.

• Non-Contact (overlap <= 0):

  Forces are typically purely normal (e.g., capillary, electrostatic), acting along \(\mathbf{n}\). The tangential force \(\mathbf{F}_t\) is zero. \(f_x = F_{\text{non-contact}} n_x\) \(f_y = F_{\text{non-contact}} n_y\) \(f_z = F_{\text{non-contact}} n_z\) (Where \(F_{\text{non-contact}}\) includes active electrostatic or capillary forces).

3. Equations Summary

• Effective Modulus:

  \[E^* = \frac{E}{2 (1 - \nu^2)}\]

• Effective Radius:

  \[R^* = \frac{R_1 R_2}{R_1 + R_2}\]

• Effective Mass:

  \[M^* = \frac{m_1 m_2}{m_1 + m_2}\]

• Overlap:

  \[\delta = (R_1 + R_2) - \text{Distance}\]

• JKR Equation (Implicit, solved for \(a\) via residual function):

  \[\delta = \frac{a^2}{R^*} - \sqrt{\frac{8 \pi a \gamma}{E^*}} \quad \Leftrightarrow \quad \text{residual}(a) = \delta - \frac{a^2}{R^*} + \sqrt{\frac{8 \pi a \gamma}{E^*}} = 0\]

• JKR Normal Force (Contact):

  \[F_n^{\text{JKR}} = \frac{4}{3} E^* \frac{a^3}{R^*} - \sqrt{6 \pi \gamma E^* a^3}\]

• Normal Damping Force:

  \[F_{n, \text{damp}} = -C_n (\mathbf{v}_{\text{rel}} \cdot \mathbf{n})\]

  where

  \[C_n = 2 \beta \sqrt{M^* k_n}, \quad \beta = \frac{-\ln(COR)}{\sqrt{\pi^2 + (\ln(COR))^2}}, \quad k_n = \frac{4}{3} E^* a\]

• Total Normal Force (Contact):

  \[F_n = F_n^{\text{JKR}} + F_{n, \text{damp}}\]

• Tangential Force (Non-History Kinetic Friction):

  \[| \mathbf{F}_t | = \mu \max(0, F_n) \quad (\text{if } |\mathbf{v}_{t, rel}| > 0 \text{ else } 0)\]

  \[\mathbf{F}_t = - | \mathbf{F}_t | \mathbf{t} \quad (\text{where } \mathbf{t} = \mathbf{v}_{t, rel} / |\mathbf{v}_{t, rel}|)\]

• Capillary Force (Non-contact, if active, as coded):

  \[F_{\text{capillary}} = 2 \pi R^* \sigma \cos(\theta) \left(1 + \min(\text{limit}, \frac{r_c}{s})\right) \quad \text{for } 0 < s < \text{s}_{Critical}\]

• Electrostatic Force (Non-contact, if active):

  \[F_{\text{elec}} = -k \frac{q_1 q_2}{\text{Distance}^2}\]

4. References

1. JKR Theory:

   • Johnson, K. L., Kendall, K., & Roberts, A. D. (1971). Surface energy and the contact of elastic solids. Proceedings of the Royal Society of London A. Mathematical and Physical Sciences, 324 (1558), 301-313.

2. Analytical Approximation (in a_return):

   • Parteli, E. J., Pöschel, T., et al. (2014). Attractive particle interaction forces and packing density of fine glass powders. Scientific Reports, 4, 6227. DOI: 10.1038/srep06227

3. General DEM:

   • Cundall, P. A., & Strack, O. D. L. (1979). A discrete numerical model for granular assemblies. Geotechnique, 29 (1), 47-65.

   • Thornton, C. (Ed.). (2015). Granular dynamics, contact mechanics and particle system simulation: a DEM approach. Springer.

4. Capillary Forces:

   • Willett, C. D., Adams, M. J., Johnson, S. A., & Seville, J. P. (2000). Capillary bridges between two spherical bodies. Langmuir, 16 (24), 9396-9405. DOI: 10.1021/la000657y