diff --git a/doc/src/pair_granular.rst b/doc/src/pair_granular.rst index 4e4d96ed06..80663b7d49 100644 --- a/doc/src/pair_granular.rst +++ b/doc/src/pair_granular.rst @@ -93,7 +93,7 @@ on particle *i* due to contact with particle *j* is given by: .. 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 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:: - \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 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 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 -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 :ref:`(Derjaguin-Muller-Toporov) ` 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} -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:: @@ -167,7 +167,7 @@ following general form: \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}`. 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` 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` -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 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 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:: - \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: @@ -294,8 +296,8 @@ keyword also affects the tangential damping. The parameter literature use :math:`x_{\gamma,t} = 1` (:ref:`Marshall `, :ref:`Tsuji et al `, :ref:`Silbert et al `). The relative 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}_r = \mathbf{v}_j - \mathbf{v}_i`. +: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` . 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 @@ -314,21 +316,24 @@ form: .. 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 :math:`F_{pulloff} = 4\pi \gamma R` for *dmt*\ . The remaining tangential options all use accumulated tangential -displacement (i.e. contact history). This is discussed below in the -context of the *linear_history* option, but the same treatment of the -accumulated displacement applies to the other options as well. +displacement (i.e. contact history), except for the options +*mindlin/force* and *mindlin_rescale/force*, that use accumulated +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: .. 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 during the entire duration of the contact: @@ -356,7 +361,7 @@ work: .. 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 current time step and :math:`\mathbf{\xi}` is the corrected @@ -372,7 +377,7 @@ discussion): .. 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 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 -For *tangential mindlin*\ , the :ref:`Mindlin ` no-slip solution is used, which differs from the *linear_history* -option by an additional factor of *a*\ , the radius of the contact region. The tangential force is given by: +For *tangential mindlin*\ , the :ref:`Mindlin ` no-slip solution +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:: - \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 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 -:math:`E` by :math:`G = E/(2(1+\nu))`, where :math:`\nu` is Poisson's ratio. This -can also be achieved by specifying *NULL* for :math:`k_t`, in which case a +discussion above. To match the Mindlin solution, one should set +:math:`k_t = 8G_{eff}`, where :math:`G_{eff}` is the effective shear modulus given by: + +.. 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 :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 -*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:: - 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 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 contact area upon unloading leads to the contact being unable to support the previous tangential loading, and spurious energy is -created without the rescaling above (:ref:`Walton ` ). See also -discussion in :ref:`Thornton et al, 2013 ` , particularly -equation 18(b) of that work and associated discussion. +created without the rescaling above (:ref:`Walton ` ). + +.. 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 ` +laws and is more consistent than *mindlin_rescale*. See discussions in +:ref:`Thornton et al, 2013 `, particularly equation 18(b) of that +work and associated discussion, and :ref:`Agnolin and Roux, 2007 `, +particularly Appendix A. ---------- @@ -460,7 +529,7 @@ exceeds a critical value: .. 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 the pseudo-force. As with tangential displacement, the rolling @@ -512,7 +581,7 @@ is then truncated according to: .. 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 displacement is rescaled so that it corresponds to the critical value @@ -763,3 +832,15 @@ Technology, 233, 30-46. .. _WaltonPC: **(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. diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index 2bf3e917cf..443478a9fc 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -43,6 +43,7 @@ Afshar agilio Agilio agni +Agnolin Ai Aidan aij @@ -599,6 +600,7 @@ Dequidt der dereference derekt +Deresiewicz Derjagin Derjaguin Derlet @@ -2219,6 +2221,7 @@ oxdna oxrna oxDNA oxRNA +packings padua Padua palegoldenrod diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index 80f208a1a2..8bf2d1bf17 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -55,7 +55,8 @@ using namespace MathSpecial; enum {HOOKE, HERTZ, HERTZ_MATERIAL, DMT, JKR}; enum {VELOCITY, MASS_VELOCITY, VISCOELASTIC, TSUJI}; 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 {ROLL_NONE, ROLL_SDS}; @@ -175,7 +176,8 @@ void PairGranular::compute(int eflag, int vflag) double signtwist, magtwist, magtortwist, Mtcrit; double tortwist1, tortwist2, tortwist3; - double shrmag,rsht; + double shrmag,rsht,prjmag; + bool frameupdate; int *ilist,*jlist,*numneigh,**firstneigh; int *touch,**firsttouch; double *history,*allhistory,**firsthistory; @@ -372,6 +374,12 @@ void PairGranular::compute(int eflag, int vflag) // 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 vt1 = vr1 - vn1; vt2 = vr2 - vn2; @@ -414,12 +422,15 @@ void PairGranular::compute(int eflag, int vflag) damp_normal_prefactor; 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; } else if (tangential_model[itype][jtype] == - TANGENTIAL_MINDLIN_RESCALE) { + TANGENTIAL_MINDLIN_RESCALE || + tangential_model[itype][jtype] == + TANGENTIAL_MINDLIN_RESCALE_FORCE) { k_tangential *= a; - // on unloading, rescale the shear displacements + // on unloading, rescale the shear displacements/force if (a < history[3]) { double factor = a/history[3]; history[0] *= factor; @@ -427,37 +438,66 @@ void PairGranular::compute(int eflag, int vflag) history[2] *= factor; } } - // rotate and update displacements. + // rotate and update displacements / force. // see e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235 if (historyupdate) { rsht = history[0]*nx + history[1]*ny + history[2]*nz; - if (fabs(rsht) < EPSILON) rsht = 0; - if (rsht > 0) { + if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_FORCE || + 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] + history[2]*history[2]); - // if rsht == shrmag, contacting pair has rotated 90 deg - // in one step, in which case you deserve a crash! - scalefac = shrmag/(shrmag - rsht); + // projection history[0] -= rsht*nx; history[1] -= rsht*ny; history[2] -= rsht*nz; + // 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[1] *= scalefac; history[2] *= scalefac; } // update history - history[0] += vtr1*dt; - history[1] += vtr2*dt; - history[2] += vtr3*dt; - if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) + if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY || + tangential_model[itype][jtype] == TANGENTIAL_MINDLIN || + 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; } // tangential forces = history + tangential velocity damping - fs1 = -k_tangential*history[0] - damp_tangential*vtr1; - fs2 = -k_tangential*history[1] - damp_tangential*vtr2; - fs3 = -k_tangential*history[2] - damp_tangential*vtr3; + if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY || + tangential_model[itype][jtype] == TANGENTIAL_MINDLIN || + 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 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] + history[2]*history[2]); if (shrmag != 0.0) { - history[0] = -1.0/k_tangential*(Fscrit*fs1/fs + - 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); + if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY || + tangential_model[itype][jtype] == TANGENTIAL_MINDLIN || + tangential_model[itype][jtype] == + TANGENTIAL_MINDLIN_RESCALE) { + history[0] = -1.0/k_tangential*(Fscrit*fs1/fs + + 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; fs2 *= Fscrit/fs; fs3 *= Fscrit/fs; @@ -512,18 +561,27 @@ void PairGranular::compute(int eflag, int vflag) int rhist1 = rhist0 + 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 (fabs(rolldotn) < EPSILON) rolldotn = 0; - if (rolldotn > 0) { // rotate into tangential plane + rolldotn = history[rhist0]*nx + history[rhist1]*ny + history[rhist2]*nz; + frameupdate = fabs(rolldotn)*k_roll < EPSILON*Frcrit; + if (frameupdate) { // rotate into tangential plane rollmag = sqrt(history[rhist0]*history[rhist0] + history[rhist1]*history[rhist1] + history[rhist2]*history[rhist2]); - scalefac = rollmag/(rollmag - rolldotn); + // projection history[rhist0] -= rolldotn*nx; history[rhist1] -= rolldotn*ny; history[rhist2] -= rolldotn*nz; // 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[rhist1] *= scalefac; history[rhist2] *= scalefac; @@ -533,14 +591,11 @@ void PairGranular::compute(int eflag, int vflag) 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; fr2 = -k_roll*history[rhist1] - damp_roll*vrl2; fr3 = -k_roll*history[rhist2] - damp_roll*vrl3; // rescale frictional displacements and forces if needed - Frcrit = roll_coeffs[itype][jtype][2] * Fncrit; fr = sqrt(fr1*fr1 + fr2*fr2 + fr3*fr3); if (fr > Frcrit) { @@ -734,7 +789,8 @@ void PairGranular::coeff(int narg, char **arg) //Defaults 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; int iarg = 2; @@ -820,7 +876,9 @@ void PairGranular::coeff(int narg, char **arg) iarg += 4; } else if ((strcmp(arg[iarg+1], "linear_history") == 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) error->all(FLERR,"Illegal pair_coeff command, " "not enough parameters provided for tangential model"); @@ -830,9 +888,15 @@ void PairGranular::coeff(int narg, char **arg) tangential_model_one = TANGENTIAL_MINDLIN; else if (strcmp(arg[iarg+1], "mindlin_rescale") == 0) 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; 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)) { if (normal_model_one == HERTZ || normal_model_one == HOOKE) { 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 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; roll_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 //**************************************** + // For linear, mindlin, mindlin_rescale: + // history = cumulative tangential displacement + // + // For mindlin/force, mindlin_rescale/force: + // history = cumulative tangential elastic force + // tangential component vt1 = vr1 - vn1; 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; if (tangential_history) { - if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN) { - k_tangential *= a; - } else if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) { + if (tangential_model[itype][jtype] != TANGENTIAL_HISTORY) { k_tangential *= a; } @@ -1529,9 +1598,17 @@ double PairGranular::single(int i, int j, int itype, int jtype, history[2]*history[2]); // tangential forces = history + tangential velocity damping - fs1 = -k_tangential*history[0] - damp_tangential*vtr1; - fs2 = -k_tangential*history[1] - damp_tangential*vtr2; - fs3 = -k_tangential*history[2] - damp_tangential*vtr3; + if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY || + tangential_model[itype][jtype] == TANGENTIAL_MINDLIN || + 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 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; fs2 *= Fscrit/fs; fs3 *= Fscrit/fs; - fs *= Fscrit/fs; + fs *= Fscrit/fs; } else fs1 = fs2 = fs3 = fs = 0.0; }