ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #2196 from jibril-b-coulibaly/mindlin_rescale

Browse Source 
Implement force history in Mindlin granular pair styles
This commit is contained in:
Axel Kohlmeyer
2020-08-26 11:42:06 -04:00
committed by GitHub
parent e2fd95d8ed 921b6d8135
commit 4ad68c98a0
3 changed files with 234 additions and 73 deletions
Download Patch File Download Diff File Expand all files Collapse all files

143
doc/src/pair_granular.rst
View File

@ -93,7 +93,7 @@ on particle *i* due to contact with particle *j* is given by:
.. math:: .. math::
\mathbf{F}_{ne, Hooke} = k_N \delta_{ij} \mathbf{n} \mathbf{F}_{ne, Hooke} = k_n \delta_{ij} \mathbf{n}
Where :math:`\delta_{ij} = R_i + R_j - \|\mathbf{r}_{ij}\|` is the particle Where :math:`\delta_{ij} = R_i + R_j - \|\mathbf{r}_{ij}\|` is the particle
overlap, :math:`R_i, R_j` are the particle radii, :math:`\mathbf{r}_{ij} = \mathbf{r}_i - \mathbf{r}_j` is the vector separating the two overlap, :math:`R_i, R_j` are the particle radii, :math:`\mathbf{r}_{ij} = \mathbf{r}_i - \mathbf{r}_j` is the vector separating the two
@ -106,7 +106,7 @@ For the *hertz* model, the normal component of force is given by:
.. math:: .. math::
\mathbf{F}_{ne, Hertz} = k_N R_{eff}^{1/2}\delta_{ij}^{3/2} \mathbf{n} \mathbf{F}_{ne, Hertz} = k_n R_{eff}^{1/2}\delta_{ij}^{3/2} \mathbf{n}
Here, :math:`R_{eff} = \frac{R_i R_j}{R_i + R_j}` is the effective Here, :math:`R_{eff} = \frac{R_i R_j}{R_i + R_j}` is the effective
radius, denoted for simplicity as *R* from here on. For *hertz*\ , the radius, denoted for simplicity as *R* from here on. For *hertz*\ , the
@ -123,7 +123,7 @@ Here, :math:`E_{eff} = E = \left(\frac{1-\nu_i^2}{E_i} + \frac{1-\nu_j^2}{E_j}\r
modulus, with :math:`\nu_i, \nu_j` the Poisson ratios of the particles of modulus, with :math:`\nu_i, \nu_j` the Poisson ratios of the particles of
types *i* and *j*\ . Note that if the elastic modulus and the shear types *i* and *j*\ . Note that if the elastic modulus and the shear
modulus of the two particles are the same, the *hertz/material* model modulus of the two particles are the same, the *hertz/material* model
is equivalent to the *hertz* model with :math:`k_N = 4/3 E_{eff}` is equivalent to the *hertz* model with :math:`k_n = 4/3 E_{eff}`
The *dmt* model corresponds to the The *dmt* model corresponds to the
:ref:`(Derjaguin-Muller-Toporov) <DMT1975>` cohesive model, where the force :ref:`(Derjaguin-Muller-Toporov) <DMT1975>` cohesive model, where the force
@ -140,7 +140,7 @@ where the force is computed as:
\mathbf{F}_{ne, jkr} = \left(\frac{4Ea^3}{3R} - 2\pi a^2\sqrt{\frac{4\gamma E}{\pi a}}\right)\mathbf{n} \mathbf{F}_{ne, jkr} = \left(\frac{4Ea^3}{3R} - 2\pi a^2\sqrt{\frac{4\gamma E}{\pi a}}\right)\mathbf{n}
Here, *a* is the radius of the contact zone, related to the overlap Here, :math:`a` is the radius of the contact zone, related to the overlap
:math:`\delta` according to: :math:`\delta` according to:
.. math:: .. math::
@ -167,7 +167,7 @@ following general form:
\mathbf{F}_{n,damp} = -\eta_n \mathbf{v}_{n,rel} \mathbf{F}_{n,damp} = -\eta_n \mathbf{v}_{n,rel}
Here, :math:`\mathbf{v}_{n,rel} = (\mathbf{v}_j - \mathbf{v}_i) \cdot \mathbf{n} \mathbf{n}` is the component of relative velocity along Here, :math:`\mathbf{v}_{n,rel} = (\mathbf{v}_j - \mathbf{v}_i) \cdot \mathbf{n}\ \mathbf{n}` is the component of relative velocity along
:math:`\mathbf{n}`. :math:`\mathbf{n}`.
The optional *damping* keyword to the *pair_coeff* command followed by The optional *damping* keyword to the *pair_coeff* command followed by
@ -259,7 +259,9 @@ tangential model choices and their expected parameters are as follows:
1. *linear_nohistory* : :math:`x_{\gamma,t}`, :math:`\mu_s` 1. *linear_nohistory* : :math:`x_{\gamma,t}`, :math:`\mu_s`
2. *linear_history* : :math:`k_t`, :math:`x_{\gamma,t}`, :math:`\mu_s` 2. *linear_history* : :math:`k_t`, :math:`x_{\gamma,t}`, :math:`\mu_s`
3. *mindlin* : :math:`k_t` or NULL, :math:`x_{\gamma,t}`, :math:`\mu_s` 3. *mindlin* : :math:`k_t` or NULL, :math:`x_{\gamma,t}`, :math:`\mu_s`
4. *mindlin_rescale* : :math:`k_t` or NULL, :math:`x_{\gamma,t}`, :math:`\mu_s` 4. *mindlin/force* : :math:`k_t` or NULL, :math:`x_{\gamma,t}`, :math:`\mu_s`
5. *mindlin_rescale* : :math:`k_t` or NULL, :math:`x_{\gamma,t}`, :math:`\mu_s`
6. *mindlin_rescale/force* : :math:`k_t` or NULL, :math:`x_{\gamma,t}`, :math:`\mu_s`
Here, :math:`x_{\gamma,t}` is a dimensionless multiplier for the normal Here, :math:`x_{\gamma,t}` is a dimensionless multiplier for the normal
damping :math:`\eta_n` that determines the magnitude of the tangential damping :math:`\eta_n` that determines the magnitude of the tangential
@ -268,11 +270,11 @@ coefficient, and :math:`k_t` is the tangential stiffness coefficient.
For *tangential linear_nohistory*, a simple velocity-dependent Coulomb For *tangential linear_nohistory*, a simple velocity-dependent Coulomb
friction criterion is used, which mimics the behavior of the *pair friction criterion is used, which mimics the behavior of the *pair
gran/hooke* style. The tangential force (\mathbf{F}_t\) is given by: gran/hooke* style. The tangential force :math:`\mathbf{F}_t` is given by:
.. math:: .. math::
\mathbf{F}_t = -min(\mu_t F_{n0}, \|\mathbf{F}_\mathrm{t,damp}\|) \mathbf{t} \mathbf{F}_t = -\min(\mu_t F_{n0}, \|\mathbf{F}_\mathrm{t,damp}\|) \mathbf{t}
The tangential damping force :math:`\mathbf{F}_\mathrm{t,damp}` is given by: The tangential damping force :math:`\mathbf{F}_\mathrm{t,damp}` is given by:
@ -294,8 +296,8 @@ keyword also affects the tangential damping. The parameter
literature use :math:`x_{\gamma,t} = 1` (:ref:`Marshall <Marshall2009>`, literature use :math:`x_{\gamma,t} = 1` (:ref:`Marshall <Marshall2009>`,
:ref:`Tsuji et al <Tsuji1992>`, :ref:`Silbert et al <Silbert2001>`). The relative :ref:`Tsuji et al <Tsuji1992>`, :ref:`Silbert et al <Silbert2001>`). The relative
tangential velocity at the point of contact is given by tangential velocity at the point of contact is given by
:math:`\mathbf{v}_{t, rel} = \mathbf{v}_{t} - (R_i\Omega_i + R_j\Omega_j) \times \mathbf{n}`, where :math:`\mathbf{v}_{t} = \mathbf{v}_r - \mathbf{v}_r\cdot\mathbf{n}{n}`, :math:`\mathbf{v}_{t, rel} = \mathbf{v}_{t} - (R_i\mathbf{\Omega}_i + R_j\mathbf{\Omega}_j) \times \mathbf{n}`, where :math:`\mathbf{v}_{t} = \mathbf{v}_r - \mathbf{v}_r\cdot\mathbf{n}\ \mathbf{n}`,
:math:`\mathbf{v}_r = \mathbf{v}_j - \mathbf{v}_i`. :math:`\mathbf{v}_r = \mathbf{v}_j - \mathbf{v}_i` .
The direction of the applied force is :math:`\mathbf{t} = \mathbf{v_{t,rel}}/\|\mathbf{v_{t,rel}}\|` . The direction of the applied force is :math:`\mathbf{t} = \mathbf{v_{t,rel}}/\|\mathbf{v_{t,rel}}\|` .
The normal force value :math:`F_{n0}` used to compute the critical force The normal force value :math:`F_{n0}` used to compute the critical force
@ -314,21 +316,24 @@ form:
.. math:: .. math::
F_{n0} = \|\mathbf{F}_ne + 2 F_{pulloff}\| F_{n0} = \|\mathbf{F}_{ne} + 2 F_{pulloff}\|
Where :math:`F_{pulloff} = 3\pi \gamma R` for *jkr*\ , and Where :math:`F_{pulloff} = 3\pi \gamma R` for *jkr*\ , and
:math:`F_{pulloff} = 4\pi \gamma R` for *dmt*\ . :math:`F_{pulloff} = 4\pi \gamma R` for *dmt*\ .
The remaining tangential options all use accumulated tangential The remaining tangential options all use accumulated tangential
displacement (i.e. contact history). This is discussed below in the displacement (i.e. contact history), except for the options
context of the *linear_history* option, but the same treatment of the *mindlin/force* and *mindlin_rescale/force*, that use accumulated
accumulated displacement applies to the other options as well. tangential force instead, and are discussed further below.
The accumulated tangential displacement is discussed in details below
in the context of the *linear_history* option. The same treatment of
the accumulated displacement applies to the other options as well.
For *tangential linear_history*, the tangential force is given by: For *tangential linear_history*, the tangential force is given by:
.. math:: .. math::
\mathbf{F}_t = -min(\mu_t F_{n0}, \|-k_t\mathbf{\xi} + \mathbf{F}_\mathrm{t,damp}\|) \mathbf{t} \mathbf{F}_t = -\min(\mu_t F_{n0}, \|-k_t\mathbf{\xi} + \mathbf{F}_\mathrm{t,damp}\|) \mathbf{t}
Here, :math:`\mathbf{\xi}` is the tangential displacement accumulated Here, :math:`\mathbf{\xi}` is the tangential displacement accumulated
during the entire duration of the contact: during the entire duration of the contact:
@ -356,7 +361,7 @@ work:
.. math:: .. math::
\mathbf{\xi} = \left(\mathbf{\xi'} - (\mathbf{n} \cdot \mathbf{\xi'})\mathbf{n}\right) \frac{\|\mathbf{\xi'}\|}{\|\mathbf{\xi'}\| - \mathbf{n}\cdot\mathbf{\xi'}} \mathbf{\xi} = \left(\mathbf{\xi'} - (\mathbf{n} \cdot \mathbf{\xi'})\mathbf{n}\right) \frac{\|\mathbf{\xi'}\|}{\|\mathbf{\xi'} - (\mathbf{n}\cdot\mathbf{\xi'})\mathbf{n}\|}
Here, :math:`\mathbf{\xi'}` is the accumulated displacement prior to the Here, :math:`\mathbf{\xi'}` is the accumulated displacement prior to the
current time step and :math:`\mathbf{\xi}` is the corrected current time step and :math:`\mathbf{\xi}` is the corrected
@ -372,7 +377,7 @@ discussion):
.. math:: .. math::
\mathbf{\xi} = -\frac{1}{k_t}\left(\mu_t F_{n0}\mathbf{t} + \mathbf{F}_{t,damp}\right) \mathbf{\xi} = -\frac{1}{k_t}\left(\mu_t F_{n0}\mathbf{t} - \mathbf{F}_{t,damp}\right)
The tangential force is added to the total normal force (elastic plus The tangential force is added to the total normal force (elastic plus
damping) to produce the total force on the particle. The tangential damping) to produce the total force on the particle. The tangential
@ -387,27 +392,68 @@ overlap region) to induce a torque on each particle according to:
\mathbf{\tau}_j = -(R_j - 0.5 \delta) \mathbf{n} \times \mathbf{F}_t \mathbf{\tau}_j = -(R_j - 0.5 \delta) \mathbf{n} \times \mathbf{F}_t
For *tangential mindlin*\ , the :ref:`Mindlin <Mindlin1949>` no-slip solution is used, which differs from the *linear_history* For *tangential mindlin*\ , the :ref:`Mindlin <Mindlin1949>` no-slip solution
option by an additional factor of *a*\ , the radius of the contact region. The tangential force is given by: is used which differs from the *linear_history* option by an additional factor
of :math:`a`, the radius of the contact region. The tangential force is given by:
.. math:: .. math::
\mathbf{F}_t = -min(\mu_t F_{n0}, \|-k_t a \mathbf{\xi} + \mathbf{F}_\mathrm{t,damp}\|) \mathbf{t} \mathbf{F}_t = -\min(\mu_t F_{n0}, \|-k_t a \mathbf{\xi} + \mathbf{F}_\mathrm{t,damp}\|) \mathbf{t}
Here, *a* is the radius of the contact region, given by :math:`a =\sqrt{R\delta}`
Here, :math:`a` is the radius of the contact region, given by :math:`a =\sqrt{R\delta}`
for all normal contact models, except for *jkr*\ , where it is given for all normal contact models, except for *jkr*\ , where it is given
implicitly by :math:`\delta = a^2/R - 2\sqrt{\pi \gamma a/E}`, see implicitly by :math:`\delta = a^2/R - 2\sqrt{\pi \gamma a/E}`, see
discussion above. To match the Mindlin solution, one should set :math:`k_t = 4G/(2-\nu)`, where :math:`G` is the shear modulus, related to Young's modulus discussion above. To match the Mindlin solution, one should set
:math:`E` by :math:`G = E/(2(1+\nu))`, where :math:`\nu` is Poisson's ratio. This :math:`k_t = 8G_{eff}`, where :math:`G_{eff}` is the effective shear modulus given by:
can also be achieved by specifying *NULL* for :math:`k_t`, in which case a
.. math::
G_{eff} = \left(\frac{2-\nu_i}{G_i} + \frac{2-\nu_j}{G_j}\right)^{-1}
where :math:`G` is the shear modulus, related to Young's modulus :math:`E`
and Poisson's ratio :math:`\nu` by :math:`G = E/(2(1+\nu))`. This can also be
achieved by specifying *NULL* for :math:`k_t`, in which case a
normal contact model that specifies material parameters :math:`E` and normal contact model that specifies material parameters :math:`E` and
:math:`\nu` is required (e.g. *hertz/material*\ , *dmt* or *jkr*\ ). In this :math:`\nu` is required (e.g. *hertz/material*\ , *dmt* or *jkr*\ ). In this
case, mixing of the shear modulus for different particle types *i* and case, mixing of the shear modulus for different particle types *i* and
*j* is done according to: *j* is done according to the formula above.
.. note::
The radius of the contact region :math:`a` depends on the normal overlap.
As a result, the tangential force for *mindlin* can change due to
a variation in normal overlap, even with no change in tangential displacement.
For *tangential mindlin/force*, the accumulated elastic tangential force
characterizes the contact history, instead of the accumulated tangential
displacement. This prevents the dependence of the tangential force on the
normal overlap as noted above. The tangential force is given by:
.. math:: .. math::
1/G = 2(2-\nu_i)(1+\nu_i)/E_i + 2(2-\nu_j)(1+\nu_j)/E_j \mathbf{F}_t = -\min(\mu_t F_{n0}, \|\mathbf{F}_{te} + \mathbf{F}_\mathrm{t,damp}\|) \mathbf{t}
The increment of the elastic component of the tangential force
:math:`\mathbf{F}_{te}` is given by:
.. math::
\mathrm{d}\mathbf{F}_{te} = -k_t a \mathbf{v}_{t,rel} \mathrm{d}\tau
The changes in frame of reference of the contacting pair of particles during
contact are accounted for by the same formula as above, replacing the
accumulated tangential displacement :math:`\xi`, by the accumulated tangential
elastic force :math:`F_{te}`. When the tangential force exceeds the critical
force, the tangential force is directly re-scaled to match the value for
the critical force:
.. math::
\mathbf{F}_{te} = - \mu_t F_{n0}\mathbf{t} + \mathbf{F}_{t,damp}
The same rules as those described for *mindlin* apply regarding the tangential
stiffness and mixing of the shear modulus for different particle types.
The *mindlin_rescale* option uses the same form as *mindlin*\ , but the The *mindlin_rescale* option uses the same form as *mindlin*\ , but the
magnitude of the tangential displacement is re-scaled as the contact magnitude of the tangential displacement is re-scaled as the contact
@ -421,9 +467,32 @@ Here, :math:`t_{n-1}` indicates the value at the previous time
step. This rescaling accounts for the fact that a decrease in the step. This rescaling accounts for the fact that a decrease in the
contact area upon unloading leads to the contact being unable to contact area upon unloading leads to the contact being unable to
support the previous tangential loading, and spurious energy is support the previous tangential loading, and spurious energy is
created without the rescaling above (:ref:`Walton <WaltonPC>` ). See also created without the rescaling above (:ref:`Walton <WaltonPC>` ).
discussion in :ref:`Thornton et al, 2013 <Thornton2013>` , particularly
equation 18(b) of that work and associated discussion. .. note::
For *mindlin*, a decrease in the tangential force already occurs as the
contact unloads, due to the dependence of the tangential force on the normal
force described above. By re-scaling :math:`\xi`, *mindlin_rescale*
effectively re-scales the tangential force twice, i.e., proportionally to
:math:`a^2`. This peculiar behavior results from use of the accumulated
tangential displacement to characterize the contact history. Although
*mindlin_rescale* remains available for historic reasons and backward
compatibility purposes, it should be avoided in favor of *mindlin_rescale/force*.
The *mindlin_rescale/force* option uses the same form as *mindlin/force*,
but the magnitude of the tangential elastic force is re-scaled as the contact
unloads, i.e. if :math:`a < a_{t_{n-1}}`:
.. math::
\mathbf{F}_{te} = \mathbf{F}_{te, t_{n-1}} \frac{a}{a_{t_{n-1}}}
This approach provides a better approximation of the :ref:`Mindlin-Deresiewicz <Mindlin1953>`
laws and is more consistent than *mindlin_rescale*. See discussions in
:ref:`Thornton et al, 2013 <Thornton2013>`, particularly equation 18(b) of that
work and associated discussion, and :ref:`Agnolin and Roux, 2007 <AgnolinRoux2007>`,
particularly Appendix A.
---------- ----------
@ -460,7 +529,7 @@ exceeds a critical value:
.. math:: .. math::
\mathbf{F}_{roll} = min(\mu_{roll} F_{n,0}, \|\mathbf{F}_{roll,0}\|)\mathbf{k} \mathbf{F}_{roll} = \min(\mu_{roll} F_{n,0}, \|\mathbf{F}_{roll,0}\|)\mathbf{k}
Here, :math:`\mathbf{k} = \mathbf{v}_{roll}/\|\mathbf{v}_{roll}\|` is the direction of Here, :math:`\mathbf{k} = \mathbf{v}_{roll}/\|\mathbf{v}_{roll}\|` is the direction of
the pseudo-force. As with tangential displacement, the rolling the pseudo-force. As with tangential displacement, the rolling
@ -512,7 +581,7 @@ is then truncated according to:
.. math:: .. math::
\tau_{twist} = min(\mu_{twist} F_{n,0}, \tau_{twist,0}) \tau_{twist} = \min(\mu_{twist} F_{n,0}, \tau_{twist,0})
Similar to the sliding and rolling displacement, the angular Similar to the sliding and rolling displacement, the angular
displacement is rescaled so that it corresponds to the critical value displacement is rescaled so that it corresponds to the critical value
@ -763,3 +832,15 @@ Technology, 233, 30-46.
.. _WaltonPC: .. _WaltonPC:
**(Otis R. Walton)** Walton, O.R., Personal Communication **(Otis R. Walton)** Walton, O.R., Personal Communication
.. _Mindlin1953:
**(Mindlin and Deresiewicz, 1953)** Mindlin, R.D., & Deresiewicz, H (1953).
Elastic Spheres in Contact under Varying Oblique Force.
J. Appl. Mech., ASME 20, 327-344.
.. _AgnolinRoux2007:
**(Agnolin and Roux 2007)** Agnolin, I. & Roux, J-N. (2007).
Internal states of model isotropic granular packings.
I. Assembling process, geometry, and contact networks. Phys. Rev. E, 76, 061302.

3
doc/utils/sphinx-config/false_positives.txt
View File

@ -43,6 +43,7 @@ Afshar
agilio agilio
Agilio Agilio
agni agni
Agnolin
Ai Ai
Aidan Aidan
aij aij
@ -599,6 +600,7 @@ Dequidt
der der
dereference dereference
derekt derekt
Deresiewicz
Derjagin Derjagin
Derjaguin Derjaguin
Derlet Derlet
@ -2219,6 +2221,7 @@ oxdna
oxrna oxrna
oxDNA oxDNA
oxRNA oxRNA
packings
padua padua
Padua Padua
palegoldenrod palegoldenrod

161
src/GRANULAR/pair_granular.cpp
View File

@ -55,7 +55,8 @@ using namespace MathSpecial;
enum {HOOKE, HERTZ, HERTZ_MATERIAL, DMT, JKR}; enum {HOOKE, HERTZ, HERTZ_MATERIAL, DMT, JKR};
enum {VELOCITY, MASS_VELOCITY, VISCOELASTIC, TSUJI}; enum {VELOCITY, MASS_VELOCITY, VISCOELASTIC, TSUJI};
enum {TANGENTIAL_NOHISTORY, TANGENTIAL_HISTORY, enum {TANGENTIAL_NOHISTORY, TANGENTIAL_HISTORY,
TANGENTIAL_MINDLIN, TANGENTIAL_MINDLIN_RESCALE}; TANGENTIAL_MINDLIN, TANGENTIAL_MINDLIN_RESCALE,
TANGENTIAL_MINDLIN_FORCE, TANGENTIAL_MINDLIN_RESCALE_FORCE};
enum {TWIST_NONE, TWIST_SDS, TWIST_MARSHALL}; enum {TWIST_NONE, TWIST_SDS, TWIST_MARSHALL};
enum {ROLL_NONE, ROLL_SDS}; enum {ROLL_NONE, ROLL_SDS};
@ -175,7 +176,8 @@ void PairGranular::compute(int eflag, int vflag)
double signtwist, magtwist, magtortwist, Mtcrit; double signtwist, magtwist, magtortwist, Mtcrit;
double tortwist1, tortwist2, tortwist3; double tortwist1, tortwist2, tortwist3;
double shrmag,rsht; double shrmag,rsht,prjmag;
bool frameupdate;
int *ilist,*jlist,*numneigh,**firstneigh; int *ilist,*jlist,*numneigh,**firstneigh;
int *touch,**firsttouch; int *touch,**firsttouch;
double *history,*allhistory,**firsthistory; double *history,*allhistory,**firsthistory;
@ -372,6 +374,12 @@ void PairGranular::compute(int eflag, int vflag)
// tangential force, including history effects // tangential force, including history effects
//**************************************** //****************************************
// For linear, mindlin, mindlin_rescale:
// history = cumulative tangential displacement
//
// For mindlin/force, mindlin_rescale/force:
// history = cumulative tangential elastic force
// tangential component // tangential component
vt1 = vr1 - vn1; vt1 = vr1 - vn1;
vt2 = vr2 - vn2; vt2 = vr2 - vn2;
@ -414,12 +422,15 @@ void PairGranular::compute(int eflag, int vflag)
damp_normal_prefactor; damp_normal_prefactor;
if (tangential_history) { if (tangential_history) {
if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN) { if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN ||
tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_FORCE) {
k_tangential *= a; k_tangential *= a;
} else if (tangential_model[itype][jtype] == } else if (tangential_model[itype][jtype] ==
TANGENTIAL_MINDLIN_RESCALE) { TANGENTIAL_MINDLIN_RESCALE ||
tangential_model[itype][jtype] ==
TANGENTIAL_MINDLIN_RESCALE_FORCE) {
k_tangential *= a; k_tangential *= a;
// on unloading, rescale the shear displacements // on unloading, rescale the shear displacements/force
if (a < history[3]) { if (a < history[3]) {
double factor = a/history[3]; double factor = a/history[3];
history[0] *= factor; history[0] *= factor;
@ -427,37 +438,66 @@ void PairGranular::compute(int eflag, int vflag)
history[2] *= factor; history[2] *= factor;
} }
} }
// rotate and update displacements. // rotate and update displacements / force.
// see e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235 // see e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235
if (historyupdate) { if (historyupdate) {
rsht = history[0]*nx + history[1]*ny + history[2]*nz; rsht = history[0]*nx + history[1]*ny + history[2]*nz;
if (fabs(rsht) < EPSILON) rsht = 0; if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_FORCE ||
if (rsht > 0) { tangential_model[itype][jtype] ==
TANGENTIAL_MINDLIN_RESCALE_FORCE)
frameupdate = fabs(rsht) < EPSILON*Fscrit;
else
frameupdate = fabs(rsht)*k_tangential < EPSILON*Fscrit;
if (frameupdate) {
shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + shrmag = sqrt(history[0]*history[0] + history[1]*history[1] +
history[2]*history[2]); history[2]*history[2]);
// if rsht == shrmag, contacting pair has rotated 90 deg // projection
// in one step, in which case you deserve a crash!
scalefac = shrmag/(shrmag - rsht);
history[0] -= rsht*nx; history[0] -= rsht*nx;
history[1] -= rsht*ny; history[1] -= rsht*ny;
history[2] -= rsht*nz; history[2] -= rsht*nz;
// also rescale to preserve magnitude // also rescale to preserve magnitude
prjmag = sqrt(history[0]*history[0] + history[1]*history[1] +
history[2]*history[2]);
if (prjmag > 0) scalefac = shrmag/prjmag;
else scalefac = 0;
history[0] *= scalefac; history[0] *= scalefac;
history[1] *= scalefac; history[1] *= scalefac;
history[2] *= scalefac; history[2] *= scalefac;
} }
// update history // update history
history[0] += vtr1*dt; if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY ||
history[1] += vtr2*dt; tangential_model[itype][jtype] == TANGENTIAL_MINDLIN ||
history[2] += vtr3*dt; tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) {
if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) // tangential displacement
history[0] += vtr1*dt;
history[1] += vtr2*dt;
history[2] += vtr3*dt;
} else {
// tangential force
// see e.g. eq. 18 of Thornton et al, Pow. Tech. 2013, v223,p30-46
history[0] -= k_tangential*vtr1*dt;
history[1] -= k_tangential*vtr2*dt;
history[2] -= k_tangential*vtr3*dt;
}
if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE ||
tangential_model[itype][jtype] ==
TANGENTIAL_MINDLIN_RESCALE_FORCE)
history[3] = a; history[3] = a;
} }
// tangential forces = history + tangential velocity damping // tangential forces = history + tangential velocity damping
fs1 = -k_tangential*history[0] - damp_tangential*vtr1; if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY ||
fs2 = -k_tangential*history[1] - damp_tangential*vtr2; tangential_model[itype][jtype] == TANGENTIAL_MINDLIN ||
fs3 = -k_tangential*history[2] - damp_tangential*vtr3; tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) {
fs1 = -k_tangential*history[0] - damp_tangential*vtr1;
fs2 = -k_tangential*history[1] - damp_tangential*vtr2;
fs3 = -k_tangential*history[2] - damp_tangential*vtr3;
} else {
fs1 = history[0] - damp_tangential*vtr1;
fs2 = history[1] - damp_tangential*vtr2;
fs3 = history[2] - damp_tangential*vtr3;
}
// rescale frictional displacements and forces if needed // rescale frictional displacements and forces if needed
fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3);
@ -465,12 +505,21 @@ void PairGranular::compute(int eflag, int vflag)
shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + shrmag = sqrt(history[0]*history[0] + history[1]*history[1] +
history[2]*history[2]); history[2]*history[2]);
if (shrmag != 0.0) { if (shrmag != 0.0) {
history[0] = -1.0/k_tangential*(Fscrit*fs1/fs + if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY ||
damp_tangential*vtr1); tangential_model[itype][jtype] == TANGENTIAL_MINDLIN ||
history[1] = -1.0/k_tangential*(Fscrit*fs2/fs + tangential_model[itype][jtype] ==
damp_tangential*vtr2); TANGENTIAL_MINDLIN_RESCALE) {
history[2] = -1.0/k_tangential*(Fscrit*fs3/fs + history[0] = -1.0/k_tangential*(Fscrit*fs1/fs +
damp_tangential*vtr3); damp_tangential*vtr1);
history[1] = -1.0/k_tangential*(Fscrit*fs2/fs +
damp_tangential*vtr2);
history[2] = -1.0/k_tangential*(Fscrit*fs3/fs +
damp_tangential*vtr3);
} else {
history[0] = Fscrit*fs1/fs + damp_tangential*vtr1;
history[1] = Fscrit*fs2/fs + damp_tangential*vtr2;
history[2] = Fscrit*fs3/fs + damp_tangential*vtr3;
}
fs1 *= Fscrit/fs; fs1 *= Fscrit/fs;
fs2 *= Fscrit/fs; fs2 *= Fscrit/fs;
fs3 *= Fscrit/fs; fs3 *= Fscrit/fs;
@ -512,18 +561,27 @@ void PairGranular::compute(int eflag, int vflag)
int rhist1 = rhist0 + 1; int rhist1 = rhist0 + 1;
int rhist2 = rhist1 + 1; int rhist2 = rhist1 + 1;
rolldotn = history[rhist0]*nx + history[rhist1]*ny + history[rhist2]*nz; k_roll = roll_coeffs[itype][jtype][0];
damp_roll = roll_coeffs[itype][jtype][1];
Frcrit = roll_coeffs[itype][jtype][2] * Fncrit;
if (historyupdate) { if (historyupdate) {
if (fabs(rolldotn) < EPSILON) rolldotn = 0; rolldotn = history[rhist0]*nx + history[rhist1]*ny + history[rhist2]*nz;
if (rolldotn > 0) { // rotate into tangential plane frameupdate = fabs(rolldotn)*k_roll < EPSILON*Frcrit;
if (frameupdate) { // rotate into tangential plane
rollmag = sqrt(history[rhist0]*history[rhist0] + rollmag = sqrt(history[rhist0]*history[rhist0] +
history[rhist1]*history[rhist1] + history[rhist1]*history[rhist1] +
history[rhist2]*history[rhist2]); history[rhist2]*history[rhist2]);
scalefac = rollmag/(rollmag - rolldotn); // projection
history[rhist0] -= rolldotn*nx; history[rhist0] -= rolldotn*nx;
history[rhist1] -= rolldotn*ny; history[rhist1] -= rolldotn*ny;
history[rhist2] -= rolldotn*nz; history[rhist2] -= rolldotn*nz;
// also rescale to preserve magnitude // also rescale to preserve magnitude
prjmag = sqrt(history[rhist0]*history[rhist0] +
history[rhist1]*history[rhist1] +
history[rhist2]*history[rhist2]);
if (prjmag > 0) scalefac = rollmag/prjmag;
else scalefac = 0;
history[rhist0] *= scalefac; history[rhist0] *= scalefac;
history[rhist1] *= scalefac; history[rhist1] *= scalefac;
history[rhist2] *= scalefac; history[rhist2] *= scalefac;
@ -533,14 +591,11 @@ void PairGranular::compute(int eflag, int vflag)
history[rhist2] += vrl3*dt; history[rhist2] += vrl3*dt;
} }
k_roll = roll_coeffs[itype][jtype][0];
damp_roll = roll_coeffs[itype][jtype][1];
fr1 = -k_roll*history[rhist0] - damp_roll*vrl1; fr1 = -k_roll*history[rhist0] - damp_roll*vrl1;
fr2 = -k_roll*history[rhist1] - damp_roll*vrl2; fr2 = -k_roll*history[rhist1] - damp_roll*vrl2;
fr3 = -k_roll*history[rhist2] - damp_roll*vrl3; fr3 = -k_roll*history[rhist2] - damp_roll*vrl3;
// rescale frictional displacements and forces if needed // rescale frictional displacements and forces if needed
Frcrit = roll_coeffs[itype][jtype][2] * Fncrit;
fr = sqrt(fr1*fr1 + fr2*fr2 + fr3*fr3); fr = sqrt(fr1*fr1 + fr2*fr2 + fr3*fr3);
if (fr > Frcrit) { if (fr > Frcrit) {
@ -734,7 +789,8 @@ void PairGranular::coeff(int narg, char **arg)
//Defaults //Defaults
normal_model_one = tangential_model_one = -1; normal_model_one = tangential_model_one = -1;
roll_model_one = twist_model_one = 0; roll_model_one = ROLL_NONE;
twist_model_one = TWIST_NONE;
damping_model_one = VISCOELASTIC; damping_model_one = VISCOELASTIC;
int iarg = 2; int iarg = 2;
@ -820,7 +876,9 @@ void PairGranular::coeff(int narg, char **arg)
iarg += 4; iarg += 4;
} else if ((strcmp(arg[iarg+1], "linear_history") == 0) || } else if ((strcmp(arg[iarg+1], "linear_history") == 0) ||
(strcmp(arg[iarg+1], "mindlin") == 0) || (strcmp(arg[iarg+1], "mindlin") == 0) ||
(strcmp(arg[iarg+1], "mindlin_rescale") == 0)) { (strcmp(arg[iarg+1], "mindlin_rescale") == 0) ||
(strcmp(arg[iarg+1], "mindlin/force") == 0) ||
(strcmp(arg[iarg+1], "mindlin_rescale/force") == 0)) {
if (iarg + 4 >= narg) if (iarg + 4 >= narg)
error->all(FLERR,"Illegal pair_coeff command, " error->all(FLERR,"Illegal pair_coeff command, "
"not enough parameters provided for tangential model"); "not enough parameters provided for tangential model");
@ -830,9 +888,15 @@ void PairGranular::coeff(int narg, char **arg)
tangential_model_one = TANGENTIAL_MINDLIN; tangential_model_one = TANGENTIAL_MINDLIN;
else if (strcmp(arg[iarg+1], "mindlin_rescale") == 0) else if (strcmp(arg[iarg+1], "mindlin_rescale") == 0)
tangential_model_one = TANGENTIAL_MINDLIN_RESCALE; tangential_model_one = TANGENTIAL_MINDLIN_RESCALE;
else if (strcmp(arg[iarg+1], "mindlin/force") == 0)
tangential_model_one = TANGENTIAL_MINDLIN_FORCE;
else if (strcmp(arg[iarg+1], "mindlin_rescale/force") == 0)
tangential_model_one = TANGENTIAL_MINDLIN_RESCALE_FORCE;
tangential_history = 1; tangential_history = 1;
if ((tangential_model_one == TANGENTIAL_MINDLIN || if ((tangential_model_one == TANGENTIAL_MINDLIN ||
tangential_model_one == TANGENTIAL_MINDLIN_RESCALE) && tangential_model_one == TANGENTIAL_MINDLIN_RESCALE ||
tangential_model_one == TANGENTIAL_MINDLIN_FORCE ||
tangential_model_one == TANGENTIAL_MINDLIN_RESCALE_FORCE) &&
(strcmp(arg[iarg+2], "NULL") == 0)) { (strcmp(arg[iarg+2], "NULL") == 0)) {
if (normal_model_one == HERTZ || normal_model_one == HOOKE) { if (normal_model_one == HERTZ || normal_model_one == HOOKE) {
error->all(FLERR, "NULL setting for Mindlin tangential " error->all(FLERR, "NULL setting for Mindlin tangential "
@ -1014,7 +1078,8 @@ void PairGranular::init_style()
} }
for (int i = 1; i <= atom->ntypes; i++) for (int i = 1; i <= atom->ntypes; i++)
for (int j = i; j <= atom->ntypes; j++) for (int j = i; j <= atom->ntypes; j++)
if (tangential_model[i][j] == TANGENTIAL_MINDLIN_RESCALE) { if (tangential_model[i][j] == TANGENTIAL_MINDLIN_RESCALE ||
tangential_model[i][j] == TANGENTIAL_MINDLIN_RESCALE_FORCE) {
size_history += 1; size_history += 1;
roll_history_index += 1; roll_history_index += 1;
twist_history_index += 1; twist_history_index += 1;
@ -1484,6 +1549,12 @@ double PairGranular::single(int i, int j, int itype, int jtype,
// tangential force, including history effects // tangential force, including history effects
//**************************************** //****************************************
// For linear, mindlin, mindlin_rescale:
// history = cumulative tangential displacement
//
// For mindlin/force, mindlin_rescale/force:
// history = cumulative tangential elastic force
// tangential component // tangential component
vt1 = vr1 - vn1; vt1 = vr1 - vn1;
vt2 = vr2 - vn2; vt2 = vr2 - vn2;
@ -1519,9 +1590,7 @@ double PairGranular::single(int i, int j, int itype, int jtype,
damp_tangential = tangential_coeffs[itype][jtype][1]*damp_normal_prefactor; damp_tangential = tangential_coeffs[itype][jtype][1]*damp_normal_prefactor;
if (tangential_history) { if (tangential_history) {
if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN) { if (tangential_model[itype][jtype] != TANGENTIAL_HISTORY) {
k_tangential *= a;
} else if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) {
k_tangential *= a; k_tangential *= a;
} }
@ -1529,9 +1598,17 @@ double PairGranular::single(int i, int j, int itype, int jtype,
history[2]*history[2]); history[2]*history[2]);
// tangential forces = history + tangential velocity damping // tangential forces = history + tangential velocity damping
fs1 = -k_tangential*history[0] - damp_tangential*vtr1; if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY ||
fs2 = -k_tangential*history[1] - damp_tangential*vtr2; tangential_model[itype][jtype] == TANGENTIAL_MINDLIN ||
fs3 = -k_tangential*history[2] - damp_tangential*vtr3; tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) {
fs1 = -k_tangential*history[0] - damp_tangential*vtr1;
fs2 = -k_tangential*history[1] - damp_tangential*vtr2;
fs3 = -k_tangential*history[2] - damp_tangential*vtr3;
} else {
fs1 = history[0] - damp_tangential*vtr1;
fs2 = history[1] - damp_tangential*vtr2;
fs3 = history[2] - damp_tangential*vtr3;
}
// rescale frictional forces if needed // rescale frictional forces if needed
fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3);
@ -1540,7 +1617,7 @@ double PairGranular::single(int i, int j, int itype, int jtype,
fs1 *= Fscrit/fs; fs1 *= Fscrit/fs;
fs2 *= Fscrit/fs; fs2 *= Fscrit/fs;
fs3 *= Fscrit/fs; fs3 *= Fscrit/fs;
fs *= Fscrit/fs; fs *= Fscrit/fs;
} else fs1 = fs2 = fs3 = fs = 0.0; } else fs1 = fs2 = fs3 = fs = 0.0;
} }