Hysteretic Spring Model

The Hysteretic Spring contact model allows plastic deformation behaviors to be included in the contact mechanics equations, resulting in particles behaving in an elastic manner up to a predefined stress.

Once this stress is exceeded, the particles behave as though undergoing plastic deformation. The result is that large overlaps are achievable without excessive forces acting upon them, thus representing a compressible material. Hysteretic Spring normal force calculation is based on the Walton-Braun theory described in the works of Walton and Braun (Walton and Braun 1986), (Walton and Braun 1986).

The following figure describes the Hysteretic Spring contact model force displacement relationship.

Note: The unloading force goes to zero before the displacement recovers to the initial contacting point.

The quantity δ0 represents the residual ‘overlap’ that remains because of the plastic deformation - flattening in the contact region. The exact location of each old contact spot is not ‘remembered,’ so that particles become like new undeformed spheres once they separate. Subsequent contacts will load along a new loading path with slope, K1. Any reloading before separation follows the slope K2 until the original loading curve K1 is reached, whereupon the lower slope is followed until unloading occurs.

The normal force, FN, is defined as:

F n =   −          K 1 δ n              for  loading                           ( K 1 δ N < K 2 ( δ n − δ 0 ))   K 2 ( δ n − δ 0 )     for  unloading/reloading      ( δ n > δ 0 )           0                 for   unloading                               ( δ n ≤ δ 0 ) MathType@MTEF@5@5@+= feaahyart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaamOramaaBa aaleaacaWGUbaabeaakiabg2da9iaaykW7caaMc8UaaGPaVlabgkHi TmaacmaaeaqabeaacaaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaayk W7caaMc8UaaGPaVlaadUeadaWgaaWcbaGaaGymaaqabaGccqaH0oaz daWgaaWcbaGaamOBaaqabaGccaaMc8UaaGPaVlaaykW7caaMc8UaaG PaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaM c8UaamOzaiaad+gacaWGYbGaaGPaVlaaykW7caWGSbGaam4Baiaadg gacaWGKbGaamyAaiaad6gacaWGNbGaaGPaVlaaykW7caaMc8UaaGPa VlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8 UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7 caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8Uaaiikai aadUeadaWgaaWcbaGaaGymaaqabaGccqaH0oazdaWgaaWcbaGaamOt aaqabaGccqGH8aapcaWGlbWaaSbaaSqaaiaaikdaaeqaaOGaaiikai abes7aKnaaBaaaleaacaWGUbaabeaakiabgkHiTiabes7aKnaaBaaa leaacaaIWaaabeaakiaacMcacaGGPaaabaGaaGPaVlaadUeadaWgaa WcbaGaaGOmaaqabaGccaGGOaGaeqiTdq2aaSbaaSqaaiaad6gaaeqa aOGaeyOeI0IaeqiTdq2aaSbaaSqaaiaaicdaaeqaaOGaaiykaiaayk W7caaMc8UaaGPaVlaaykW7caaMc8UaamOzaiaad+gacaWGYbGaaGPa VlaaykW7caWG1bGaamOBaiaadYgacaWGVbGaamyyaiaadsgacaWGPb GaamOBaiaadEgacaGGVaGaamOCaiaadwgacaWGSbGaam4Baiaadgga caWGKbGaamyAaiaad6gacaWGNbGaaGPaVlaaykW7caaMc8UaaGPaVl aaykW7caaMc8Uaaiikaiabes7aKnaaBaaaleaacaWGUbaabeaakiab g6da+iabes7aKnaaBaaaleaacaaIWaaabeaakiaacMcaaeaacaaMc8 UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7 caaMc8UaaGimaiaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVl aaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8Ua aGPaVlaaykW7caaMc8UaamOzaiaad+gacaWGYbGaaGPaVlaaykW7ca aMc8UaamyDaiaad6gacaWGSbGaam4BaiaadggacaWGKbGaamyAaiaa d6gacaWGNbGaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaG PaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaM c8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaayk W7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPa VlaacIcacqaH0oazdaWgaaWcbaGaamOBaaqabaGccqGHKjYOcqaH0o azdaWgaaWcbaGaaGimaaqabaGccaGGPaaaaiaawUhacaGL9baaaaa@59A9@

Where K1 and K2 are the loading and unloading stiffness respectively, δn, is the normal overlap and δ0 is the residual overlap. Loading stiffness, K1, is related to the yield strength of each material participating in contact, Y1 and Y2, as follows (Walton O., 2006) in: 

K 1 = 5 R * min ( Y 1 , Y 2 ) MathType@MTEF@5@5@+= feaahyart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaam4samaaBa aaleaacaaIXaaabeaakiabg2da9iaaiwdacaWGsbWaaWbaaSqabeaa caGGQaaaaOGaciyBaiaacMgacaGGUbGaaiikaiaadMfadaWgaaWcba GaaGymaaqabaGccaGGSaGaamywamaaBaaaleaacaaIYaaabeaakiaa cMcaaaa@43B1@

Here, R* is the equivalent radius of the two particles participating in the contact. The Coefficient of Restitution is defined as: 

e = K 1 K 2 MathType@MTEF@5@5@+= feaahyart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaamyzaiabg2 da9maakaaabaWaaSaaaeaacaWGlbWaaSbaaSqaaiaaigdaaeqaaaGc baGaam4samaaBaaaleaacaaIYaaabeaaaaaabeaaaaa@3B7E@

This allows for the determination of K>2 from K1 (Walton O. , 2006). Residual overlap is updated every Time Step according to the following rule: 

δ n  =   −   δ n   1− K 1 K 2             for  loading                  ( K 1 δ N < K 2 ( δ n − δ 0 ))              δ 0                    for  unloading/reloading      ( δ n > δ 0 )              δ n                   for   unloading                               ( δ n ≤ δ 0 ) MathType@MTEF@5@5@+= feaahGart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaeqiTdq2aaS baaSqaaiaad6gaaeqaaOGaaGPaVlabg2da9iaaykW7caaMc8UaaGPa VlabgkHiTmaacmaaeaqabeaacaaMc8UaeqiTdq2aaSbaaSqaaiaad6 gaaeqaaOGaaGPaVpaabmaabaGaaGymaiabgkHiTmaalaaabaGaam4s amaaBaaaleaacaaIXaaabeaaaOqaaiaadUeadaWgaaWcbaGaaGOmaa qabaaaaaGccaGLOaGaayzkaaGaaGPaVlaaykW7caaMc8UaaGPaVlaa ykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8Uaam Ozaiaad+gacaWGYbGaaGPaVlaaykW7caWGSbGaam4BaiaadggacaWG KbGaamyAaiaad6gacaWGNbGaaGPaVlaaykW7caaMc8UaaGPaVlaayk W7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPa VlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaiikaiaadUeadaWgaa WcbaGaaGymaaqabaGccqaH0oazdaWgaaWcbaGaamOtaaqabaGccqGH 8aapcaWGlbWaaSbaaSqaaiaaikdaaeqaaOGaaiikaiabes7aKnaaBa aaleaacaWGUbaabeaakiabgkHiTiabes7aKnaaBaaaleaacaaIWaaa beaakiaacMcacaGGPaaabaGaaGPaVlaaykW7caaMc8UaaGPaVlaayk W7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaeqiT dq2aaSbaaSqaaiaaicdaaeqaaOGaaGPaVlaaykW7caaMc8UaaGPaVl aaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8Ua aGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaadAgaca WGVbGaamOCaiaaykW7caaMc8UaamyDaiaad6gacaWGSbGaam4Baiaa dggacaWGKbGaamyAaiaad6gacaWGNbGaai4laiaadkhacaWGLbGaam iBaiaad+gacaWGHbGaamizaiaadMgacaWGUbGaam4zaiaaykW7caaM c8UaaGPaVlaaykW7caaMc8UaaGPaVlaacIcacqaH0oazdaWgaaWcba GaamOBaaqabaGccqGH+aGpcqaH0oazdaWgaaWcbaGaaGimaaqabaGc caGGPaaabaGaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaG PaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaeqiTdq2aaSbaaSqa aiaad6gaaeqaaOGaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8 UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7 caaMc8UaaGPaVlaaykW7caaMc8UaamOzaiaad+gacaWGYbGaaGPaVl aaykW7caaMc8UaamyDaiaad6gacaWGSbGaam4BaiaadggacaWGKbGa amyAaiaad6gacaWGNbGaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7ca aMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaa ykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaG PaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaM c8UaaGPaVlaacIcacqaH0oazdaWgaaWcbaGaamOBaaqabaGccqGHKj YOcqaH0oazdaWgaaWcbaGaaGimaaqabaGccaGGPaaaaiaawUhacaGL 9baaaaa@6ED5@

The main energy dissipation mechanism is due to the difference in spring stiffness between loading and unloading phases. The need for an additional velocity dependent dissipation mechanism, to suppress possible undamped low amplitude oscillations, was suggested in (Walton O. , 2006). In EDEM's Hysteretic Spring model, this mechanism was implemented in a way similar to the linear spring normal damping force scaled by a dimensionless user set parameter Damping Factor, bn, as follows:

F n d =−b n 4 m * k 1+ π lne 2 v n rel → MathType@MTEF@5@5@+= feaahGart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaamOramaaBa aaleaacaWGUbaabeaakmaaCaaaleqabaGaamizaaaakiabg2da9iab gkHiTiaadkgacaWLa8+aaSbaaSqaaiaad6gaaeqaaOWaaOaaaeaada WcaaqaaiaaisdacaWGTbWaaWbaaSqabeaacaGGQaaaaOGaam4Aaaqa aiaaigdacqGHRaWkdaqadaqaamaalaaabaGaeqiWdahabaGaciiBai aac6gacaWGLbaaaaGaayjkaiaawMcaamaaCaaaleqabaGaaGOmaaaa aaaabeaakiaadAhadaWhcaqaamaaDaaaleaacaWGUbaabaGaamOCai aadwgacaWGSbaaaaGccaGLxdcaaaa@519A@

Here k is either K1 or K2. The rest of the variables and factors in the above expression are explained in the Linear Spring contact model section.

Since the Hysteretic Spring model for normal force can be considered as an extension of the Linear Spring model for normal force, it is recommended to use tangential components of the Linear Spring model for tangential forces. In particular, a more general form of the Linear Spring Elastic tangential force is used in the Hysteretic Spring model in which the force depends on the user-defined stiffness factor, γt, as follows:

F t = − min ( γ 1 K 1 δ t + F t d , μ F n ) MathType@MTEF@5@5@+= feaahyart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaamOramaaBa aaleaacaWG0baabeaakiabg2da9iabgkHiTiGac2gacaGGPbGaaiOB aiaacIcacqaHZoWzdaWgaaWcbaGaaGymaaqabaGccaWGlbWaaSbaaS qaaiaaigdaaeqaaOGaeqiTdq2aaSbaaSqaaiaadshaaeqaaOGaey4k aSIaamOramaaBaaaleaacaWG0baabeaakmaaCaaaleqabaGaamizaa aakiaacYcacqaH8oqBcaWGgbWaaSbaaSqaaiaad6gaaeqaaOGaaiyk aaaa@4D90@

As opposed to the normal component, the tangential damping force is not scaled by the damping factor.

F t d =− 4 m * γ t k 1+ π lne 2 v t rel → MathType@MTEF@5@5@+= feaahyart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaamOramaaBa aaleaacaWG0baabeaakmaaCaaaleqabaGaamizaaaakiabg2da9iab gkHiTmaakaaabaWaaSaaaeaacaaI0aGaamyBamaaCaaaleqabaGaai Okaaaakiabeo7aNnaaBaaaleaacaWG0baabeaakiaadUgaaeaacaaI XaGaey4kaSYaaeWaaeaadaWcaaqaaiabec8aWbqaaiGacYgacaGGUb GaamyzaaaaaiaawIcacaGLPaaadaahaaWcbeqaaiaaikdaaaaaaaqa baGccaWG2bWaa8HaaeaadaqhaaWcbaGaamiDaaqaaiaadkhacaWGLb GaamiBaaaaaOGaay51Gaaaaa@50E2@

For particles, a default value for the Yield Strength is estimated from the Young’s Modulus E and Radius of the particle R according to the algorithm suggested in (Walton O. , 2006). Initially, the possible default value is calculated as follows: 

Y =   4 15 E R 1 70 MathType@MTEF@5@5@+= feaahyart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaamywaiabg2 da9iaaykW7daWcaaqaaiaaisdaaeaacaaIXaGaaGynaaaadaWcaaqa aiaadweaaeaadaGcaaqaaiaadkfaaSqabaaaaOWaaSaaaeaacaaIXa aabaGaaG4naiaaicdaaaaaaa@3FC8@

If this value is less than 100Pa, it is replaced with a value of 0.003 *E.

For Geometries, the default value is set to 1e+09 Pa.

The following table describes the interaction and contact model parameters for the Hysteretic Spring model.
Interaction Configurable Parameters Position

Particle-Particle,

Particle-Geometry

Click the + icon to add Particle-Particle or Particle-Geometry interactions.

Set the following interaction parameters:

Damping factor

This dimensionless parameter controls the amount of velocity dependent damping. Without this additional damping, small vibrations of particles can persist for a long time. Setting this parameter to zero suppresses only the velocity-dependent damping without affecting the main energy loss mechanism due to the nature of the linearized elasto-plastic model. The recommended value for this parameter is zero when modeling “aggressive” compression of materials and small non-zero value for small stain compression regimes.

Stiffness factor

This dimensionless parameter is defined as the ratio of tangential stiffness to normal loading stiffness and used in calculation of tangential forces at the contact. According to the literature, the typical value of this parameter is in the range between 0.7 and 1.

Also set the yield strength in units of Pa for each particle or Geometry that has an interaction. EDEM offers a reasonable default value for this parameter estimated from material Shear Modulus. You can, however, overwrite this default value.

 Last