From a6cd4a935ebda9328a5c8fa2d134423f5885d707 Mon Sep 17 00:00:00 2001 From: "Jibril B. Coulibaly" Date: Sat, 27 Jun 2020 01:17:46 -0500 Subject: [PATCH 1/7] Replace accumulated displacement by accumulated force for tangential force in styles mindlin and mindlin_rescale. Change documentation accordingly --- doc/src/pair_granular.rst | 65 ++++++++++++++++++++------------ src/GRANULAR/pair_granular.cpp | 69 +++++++++++++++++++++++----------- 2 files changed, 89 insertions(+), 45 deletions(-) diff --git a/doc/src/pair_granular.rst b/doc/src/pair_granular.rst index ee90d1537a..a8509a6019 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 @@ -268,7 +268,7 @@ 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:: @@ -294,8 +294,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 @@ -319,10 +319,11 @@ form: 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. +The remaining tangential options all use accumulated tangential displacement, +or accumulated tangential force (i.e. contact history). This is discussed +in details below for accumulated tangential displacement in the context +of the *linear_history* option. The same treatment of the accumulated +displacement, or accumulated force, applies to the other options as well. For *tangential linear_history*, the tangential force is given by: @@ -372,7 +373,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,35 +388,51 @@ 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 *a*\ , the radius of the contact region. The tangential stiffness depends +on the radius of the contact region and the force must therefore be computed +incrementally. The accumulated tangential force characterizes the contact history. +The increment of the elastic component of the tangential force :math:`\mathbf{F}_{te}` +is given by: .. math:: - \mathbf{F}_t = -min(\mu_t F_{n0}, \|-k_t a \mathbf{\xi} + \mathbf{F}_\mathrm{t,damp}\|) \mathbf{t} + \mathrm{d}\mathbf{F}_{te} = -k_t a \mathbf{v}_{t,rel} \mathrm{d}\tau + +The tangential force is given by: + +.. math:: + + \mathbf{F}_t = -min(\mu_t F_{n0}, \|\mathbf{F}_{te} + \mathbf{F}_\mathrm{t,damp}\|) \mathbf{t} Here, *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 -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: +discussion above. To match the Mindlin solution, one should set +:math:`k_t = 8G_{eff}`, where :math:`G` is the effective shear modulus given by: .. math:: - 1/G = 2(2-\nu_i)(1+\nu_i)/E_i + 2(2-\nu_j)(1+\nu_j)/E_j + 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 the formula above. + + 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 force is re-scaled as the contact unloads, i.e. if :math:`a < a_{t_{n-1}}`: .. math:: - \mathbf{\xi} = \mathbf{\xi_{t_{n-1}}} \frac{a}{a_{t_{n-1}}} + \mathbf{F}_{te} = \mathbf{F}_{te, t_{n-1}} \frac{a}{a_{t_{n-1}}} Here, :math:`t_{n-1}` indicates the value at the previous time step. This rescaling accounts for the fact that a decrease in the diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index fef8ded5f7..ece1e9e211 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -377,6 +377,9 @@ void PairGranular::compute(int eflag, int vflag) // tangential force, including history effects //**************************************** + // For linear forces, history = cumulative tangential displacement + // For nonlinear forces, history = cumulative elastic tangential force + // tangential component vt1 = vr1 - vn1; vt2 = vr2 - vn2; @@ -424,7 +427,7 @@ void PairGranular::compute(int eflag, int vflag) } else if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) { k_tangential *= a; - // on unloading, rescale the shear displacements + // on unloading, rescale the shear force if (a < history[3]) { double factor = a/history[3]; history[0] *= factor; @@ -432,11 +435,11 @@ 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 (fabs(rsht) < EPSILON) rsht = 0; // EPSILON still valid when history is a force??? if (rsht > 0) { shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + history[2]*history[2]); @@ -451,18 +454,30 @@ void PairGranular::compute(int eflag, int vflag) 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) - history[3] = a; + // update tangential displacement / force history + if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY) { + history[0] += vtr1*dt; + history[1] += vtr2*dt; + history[2] += vtr3*dt; + } else { // Eq. (18) Thornton et al., 2013 (http://dx.doi.org/10.1016/j.powtec.2012.08.012) + 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) + 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) { + 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); @@ -470,12 +485,18 @@ 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) { + 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; @@ -1534,9 +1555,15 @@ 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) { + 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); From bf724332d480d80578b1fd76e614c077a86fb831 Mon Sep 17 00:00:00 2001 From: "Jibril B. Coulibaly" Date: Mon, 10 Aug 2020 10:53:30 -0500 Subject: [PATCH 2/7] implement tangential force history in mindlin/force and mindlin_rescale/force --- doc/src/pair_granular.rst | 134 ++++++++++++++++++++++++--------- src/GRANULAR/pair_granular.cpp | 83 ++++++++++++++------ 2 files changed, 158 insertions(+), 59 deletions(-) diff --git a/doc/src/pair_granular.rst b/doc/src/pair_granular.rst index a8509a6019..ab735bdfde 100644 --- a/doc/src/pair_granular.rst +++ b/doc/src/pair_granular.rst @@ -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` +5. *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 @@ -272,7 +274,7 @@ 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,7 +296,7 @@ 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\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}_{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}}\|` . @@ -314,22 +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, -or accumulated tangential force (i.e. contact history). This is discussed -in details below for accumulated tangential displacement in the context -of the *linear_history* option. The same treatment of the accumulated -displacement, or accumulated force, applies to the other options as well. +The remaining tangential options all use accumulated tangential +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: @@ -390,27 +394,18 @@ overlap region) to induce a torque on each particle according to: 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 stiffness depends -on the radius of the contact region and the force must therefore be computed -incrementally. The accumulated tangential force characterizes the contact history. -The increment of the elastic component of the tangential force :math:`\mathbf{F}_{te}` -is given by: +of :math:`a`, the radius of the contact region. The tangential force is given by: .. math:: - \mathrm{d}\mathbf{F}_{te} = -k_t a \mathbf{v}_{t,rel} \mathrm{d}\tau - -The tangential force is given by: + \mathbf{F}_t = -\min(\mu_t F_{n0}, \|-k_t a \mathbf{\xi} + \mathbf{F}_\mathrm{t,damp}\|) \mathbf{t} -.. math:: - \mathbf{F}_t = -min(\mu_t F_{n0}, \|\mathbf{F}_{te} + \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 = 8G_{eff}`, where :math:`G` is the effective shear modulus given by: +:math:`k_t = 8G_{eff}`, where :math:`G_{eff}` is the effective shear modulus given by: .. math:: @@ -424,23 +419,80 @@ normal contact model that specifies material parameters :math:`E` and case, mixing of the shear modulus for different particle types *i* and *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:: + + \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 force is re-scaled as the contact +magnitude of the tangential displacement is re-scaled as the contact +unloads, i.e. if :math:`a < a_{t_{n-1}}`: + +.. math:: + + \mathbf{\xi} = \mathbf{\xi_{t_{n-1}}} \frac{a}{a_{t_{n-1}}} + +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 ` ). + +.. 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}}} -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. +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. ---------- @@ -477,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 @@ -529,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 @@ -794,3 +846,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. \ No newline at end of file diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index ece1e9e211..e732453757 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -54,7 +54,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}; @@ -377,8 +378,11 @@ void PairGranular::compute(int eflag, int vflag) // tangential force, including history effects //**************************************** - // For linear forces, history = cumulative tangential displacement - // For nonlinear forces, history = cumulative elastic tangential force + // 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; @@ -422,12 +426,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 force + // on unloading, rescale the shear displacements/force if (a < history[3]) { double factor = a/history[3]; history[0] *= factor; @@ -435,11 +442,11 @@ void PairGranular::compute(int eflag, int vflag) history[2] *= factor; } } - // rotate and update displacements / force + // 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; // EPSILON still valid when history is a force??? + if (fabs(rsht) < EPSILON) rsht = 0; // EPSILON = 1e-10 can fail for small granular media: if (rsht > 0) { shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + history[2]*history[2]); @@ -454,22 +461,31 @@ void PairGranular::compute(int eflag, int vflag) history[1] *= scalefac; history[2] *= scalefac; } - // update tangential displacement / force history - if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY) { + // update history + 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 { // Eq. (18) Thornton et al., 2013 (http://dx.doi.org/10.1016/j.powtec.2012.08.012) + } 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) - history[3] = a; } + 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 - if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY) { + 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; @@ -485,7 +501,10 @@ 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) { - if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY) { + 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 + @@ -760,7 +779,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; @@ -846,7 +866,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"); @@ -856,9 +878,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 " @@ -1040,7 +1068,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; @@ -1510,6 +1539,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; @@ -1545,9 +1580,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; } @@ -1555,7 +1588,9 @@ double PairGranular::single(int i, int j, int itype, int jtype, history[2]*history[2]); // tangential forces = history + tangential velocity damping - if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY) { + 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; @@ -1572,7 +1607,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; } From 2de98999c139c0694a263fb90556875add79e7e1 Mon Sep 17 00:00:00 2001 From: "Jibril B. Coulibaly" Date: Mon, 10 Aug 2020 14:51:00 -0500 Subject: [PATCH 3/7] bug fix formula for frame of reference rotation for granular tangential history --- doc/src/pair_granular.rst | 2 +- src/GRANULAR/pair_granular.cpp | 57 ++++++++++++++++++---------------- 2 files changed, 31 insertions(+), 28 deletions(-) diff --git a/doc/src/pair_granular.rst b/doc/src/pair_granular.rst index ab735bdfde..0b6e12ade7 100644 --- a/doc/src/pair_granular.rst +++ b/doc/src/pair_granular.rst @@ -361,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 diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index e732453757..9bde81b981 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -446,20 +446,20 @@ void PairGranular::compute(int eflag, int vflag) // 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; // EPSILON = 1e-10 can fail for small granular media: - if (rsht > 0) { - 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); - history[0] -= rsht*nx; - history[1] -= rsht*ny; - history[2] -= rsht*nz; - // also rescale to preserve magnitude - history[0] *= scalefac; - history[1] *= scalefac; - history[2] *= scalefac; + shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + + history[2]*history[2]); + // 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]); + scalefac = shrmag/prjmag; + history[0] *= scalefac; + history[1] *= scalefac; + history[2] *= scalefac; } // update history if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY || @@ -559,19 +559,22 @@ void PairGranular::compute(int eflag, int vflag) rolldotn = history[rhist0]*nx + history[rhist1]*ny + history[rhist2]*nz; if (historyupdate) { - if (fabs(rolldotn) < EPSILON) rolldotn = 0; - if (rolldotn > 0) { // rotate into tangential plane - rollmag = sqrt(history[rhist0]*history[rhist0] + - history[rhist1]*history[rhist1] + - history[rhist2]*history[rhist2]); - scalefac = rollmag/(rollmag - rolldotn); - history[rhist0] -= rolldotn*nx; - history[rhist1] -= rolldotn*ny; - history[rhist2] -= rolldotn*nz; - // also rescale to preserve magnitude - history[rhist0] *= scalefac; - history[rhist1] *= scalefac; - history[rhist2] *= scalefac; + // rotate into tangential plane + rollmag = sqrt(history[rhist0]*history[rhist0] + + history[rhist1]*history[rhist1] + + history[rhist2]*history[rhist2]); + // 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]); + scalefac = rollmag/prjmag; + history[rhist0] *= scalefac; + history[rhist1] *= scalefac; + history[rhist2] *= scalefac; } history[rhist0] += vrl1*dt; history[rhist1] += vrl2*dt; From 5ebac27fd5d45492b02ae1d46d1fe8f345593a9e Mon Sep 17 00:00:00 2001 From: "Jibril B. Coulibaly" Date: Mon, 10 Aug 2020 15:15:47 -0500 Subject: [PATCH 4/7] safety for division by zero in scaling of the projection --- src/GRANULAR/pair_granular.cpp | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index 9bde81b981..f59c00e90a 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -181,7 +181,7 @@ void PairGranular::compute(int eflag, int vflag) double signtwist, magtwist, magtortwist, Mtcrit; double tortwist1, tortwist2, tortwist3; - double shrmag,rsht; + double shrmag,rsht,prjmag; int *ilist,*jlist,*numneigh,**firstneigh; int *touch,**firsttouch; double *history,*allhistory,**firsthistory; @@ -456,11 +456,12 @@ void PairGranular::compute(int eflag, int vflag) // also rescale to preserve magnitude prjmag = sqrt(history[0]*history[0] + history[1]*history[1] + history[2]*history[2]); - scalefac = shrmag/prjmag; + if (prjmag > 0) scalefac = shrmag/prjmag; + else scalefac = 0; history[0] *= scalefac; history[1] *= scalefac; history[2] *= scalefac; - } + // update history if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY || tangential_model[itype][jtype] == TANGENTIAL_MINDLIN || @@ -571,11 +572,12 @@ void PairGranular::compute(int eflag, int vflag) prjmag = sqrt(history[rhist0]*history[rhist0] + history[rhist1]*history[rhist1] + history[rhist2]*history[rhist2]); - scalefac = rollmag/prjmag; + if (prjmag > 0) scalefac = rollmag/prjmag; + else scalefac = 0; history[rhist0] *= scalefac; history[rhist1] *= scalefac; history[rhist2] *= scalefac; - } + history[rhist0] += vrl1*dt; history[rhist1] += vrl2*dt; history[rhist2] += vrl3*dt; From e742ae747532dd7e12558974216dab5fd45cc4d2 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 21 Aug 2020 00:31:55 -0400 Subject: [PATCH 5/7] fix RST syntax and spelling issues in granular pair style docs --- doc/src/pair_granular.rst | 8 ++++---- doc/utils/sphinx-config/false_positives.txt | 3 +++ 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/doc/src/pair_granular.rst b/doc/src/pair_granular.rst index 40f3888e8f..d846e79e35 100644 --- a/doc/src/pair_granular.rst +++ b/doc/src/pair_granular.rst @@ -261,7 +261,7 @@ tangential model choices and their expected parameters are as follows: 3. *mindlin* : :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` -5. *mindlin_rescale/force* : :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 @@ -831,14 +831,14 @@ Technology, 233, 30-46. **(Otis R. Walton)** Walton, O.R., Personal Communication -** _Mindlin1953: +.. _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: +.. _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. \ No newline at end of file +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 From ee6ef98b9be42a4d20cce584b5acaa175666f98d Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 21 Aug 2020 00:42:57 -0400 Subject: [PATCH 6/7] remove trailing whitespace --- doc/src/pair_granular.rst | 6 +++--- src/GRANULAR/pair_granular.cpp | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/doc/src/pair_granular.rst b/doc/src/pair_granular.rst index d846e79e35..2fce2d24f1 100644 --- a/doc/src/pair_granular.rst +++ b/doc/src/pair_granular.rst @@ -410,7 +410,7 @@ discussion above. To match the Mindlin solution, one should set .. 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 @@ -476,8 +476,8 @@ created without the rescaling above (:ref:`Walton ` ). 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 + 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*, diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index b16fd40236..1b92c2686d 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -447,7 +447,7 @@ void PairGranular::compute(int eflag, int vflag) 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]); @@ -572,7 +572,7 @@ void PairGranular::compute(int eflag, int vflag) history[rhist0] *= scalefac; history[rhist1] *= scalefac; history[rhist2] *= scalefac; - + history[rhist0] += vrl1*dt; history[rhist1] += vrl2*dt; history[rhist2] += vrl3*dt; From 921b6d8135d4aa92c9050e7a33d6e05b77a94d6e Mon Sep 17 00:00:00 2001 From: "Jibril B. Coulibaly" Date: Fri, 21 Aug 2020 13:20:53 -0500 Subject: [PATCH 7/7] relative threshold for contact frame update based on tangential critical force --- src/GRANULAR/pair_granular.cpp | 84 +++++++++++++++++++--------------- 1 file changed, 47 insertions(+), 37 deletions(-) diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index 1b92c2686d..8bf2d1bf17 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -177,6 +177,7 @@ void PairGranular::compute(int eflag, int vflag) double tortwist1, tortwist2, tortwist3; double shrmag,rsht,prjmag; + bool frameupdate; int *ilist,*jlist,*numneigh,**firstneigh; int *touch,**firsttouch; double *history,*allhistory,**firsthistory; @@ -441,22 +442,29 @@ void PairGranular::compute(int eflag, int vflag) // see e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235 if (historyupdate) { rsht = history[0]*nx + history[1]*ny + history[2]*nz; - shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + - history[2]*history[2]); - // 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; + 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]); + // 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 if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY || tangential_model[itype][jtype] == TANGENTIAL_MINDLIN || @@ -553,39 +561,41 @@ 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; - if (historyupdate) { - // rotate into tangential plane - rollmag = sqrt(history[rhist0]*history[rhist0] + - history[rhist1]*history[rhist1] + - history[rhist2]*history[rhist2]); - // 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; + k_roll = roll_coeffs[itype][jtype][0]; + damp_roll = roll_coeffs[itype][jtype][1]; + Frcrit = roll_coeffs[itype][jtype][2] * Fncrit; + if (historyupdate) { + 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]); + // 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; + } history[rhist0] += vrl1*dt; history[rhist1] += vrl2*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; 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) {