Merge pull request #1631 from jibril-b-coulibaly/patch-2

Bug fixes for pair style granular

doc/src/pair_granular.txt

@@ -100,7 +100,7 @@ on particle {i} due to contact with particle {j} is given by:
 \mathbf\{F\}_\{ne, Hooke\} = k_N \delta_\{ij\} \mathbf\{n\}
 \end\{equation\}
 
-Where \(\delta = R_i + R_j - \|\mathbf\{r\}_\{ij\}\|\) is the particle
+Where \(\delta_\{ij\} = R_i + R_j - \|\mathbf\{r\}_\{ij\}\|\) is the particle
 overlap, \(R_i, R_j\) are the particle radii, \(\mathbf\{r\}_\{ij\} =
 \mathbf\{r\}_i - \mathbf\{r\}_j\) is the vector separating the two
 particle centers (note the i-j ordering so that \(F_\{ne\}\) is
@@ -177,7 +177,7 @@ following general form:
 \end\{equation\}
 
 Here, \(\mathbf\{v\}_\{n,rel\} = (\mathbf\{v\}_j - \mathbf\{v\}_i)
-\cdot \mathbf\{n\}\) is the component of relative velocity along
+\cdot \mathbf\{n\} \mathbf\{n\}\) is the component of relative velocity along
 \(\mathbf\{n\}\).
 
 The optional {damping} keyword to the {pair_coeff} command followed by
@@ -299,8 +299,8 @@ the normal damping \(\eta_n\) (see above):
 \eta_t = -x_\{\gamma,t\} \eta_n
 \end\{equation\}
 
-The normal damping prefactor \(\eta_n\) is determined by the choice of
+The normal damping prefactor \(\eta_n\) is determined by the choice
-the {damping} keyword, as discussed above.  Thus, the {damping}
+of the {damping} keyword, as discussed above.  Thus, the {damping}
 keyword also affects the tangential damping.  The parameter
 \(x_\{\gamma,t\}\) is a scaling coefficient. Several works in the
 literature use \(x_\{\gamma,t\} = 1\) ("Marshall"_#Marshall2009,
@@ -308,10 +308,10 @@ literature use \(x_\{\gamma,t\} = 1\) ("Marshall"_#Marshall2009,
 tangential velocity at the point of contact is given by
 \(\mathbf\{v\}_\{t, rel\} = \mathbf\{v\}_\{t\} - (R_i\Omega_i +
 R_j\Omega_j) \times \mathbf\{n\}\), where \(\mathbf\{v\}_\{t\} =
-\mathbf\{v\}_r - \mathbf\{v\}_r\cdot\mathbf\{n\}\), \(\mathbf\{v\}_r =
+\mathbf\{v\}_r - \mathbf\{v\}_r\cdot\mathbf\{n\}\{n\}\),
-\mathbf\{v\}_j - \mathbf\{v\}_i\). The direction of the applied force
+\(\mathbf\{v\}_r = \mathbf\{v\}_j - \mathbf\{v\}_i\).
-is \(\mathbf\{t\} =
+The direction of the applied force is \(\mathbf\{t\} =
-\mathbf\{v_\{t,rel\}\}/\|\mathbf\{v_\{t,rel\}\}\|\).
+\mathbf\{v_\{t,rel\}\}/\|\mathbf\{v_\{t,rel\}\}\|\) .
 
 The normal force value \(F_\{n0\}\) used to compute the critical force
 depends on the form of the contact model. For non-cohesive models
@@ -411,8 +411,8 @@ option by an additional factor of {a}, the radius of the contact region. The tan
 \mathbf\{F\}_t =  -min(\mu_t F_\{n0\}, \|-k_t a \mathbf\{\xi\} + \mathbf\{F\}_\mathrm\{t,damp\}\|) \mathbf\{t\}
 \end\{equation\}
 
-Here, {a} is the radius of the contact region, given by \(a = \delta
+Here, {a} is the radius of the contact region, given by \(a =\sqrt\{R\delta\}\)
-R\) 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 \(\delta = a^2/R - 2\sqrt\{\pi \gamma a/E\}\), see
 discussion above. To match the Mindlin solution, one should set \(k_t
 = 8G\), where \(G\) is the shear modulus, related to Young's modulus
@@ -680,7 +680,7 @@ The single() function of these pair styles returns 0.0 for the energy
 of a pairwise interaction, since energy is not conserved in these
 dissipative potentials.  It also returns only the normal component of
 the pairwise interaction force.  However, the single() function also
-calculates 10 extra pairwise quantities.  The first 3 are the
+calculates 12 extra pairwise quantities.  The first 3 are the
 components of the tangential force between particles I and J, acting
 on particle I.  The 4th is the magnitude of this tangential force.
 The next 3 (5-7) are the components of the rolling torque acting on

src/GRANULAR/pair_granular.cpp

@@ -305,6 +305,7 @@ void PairGranular::compute(int eflag, int vflag)
 
         delta = radsum - r;
         dR = delta*Reff;
+
 	if (normal_model[itype][jtype] == JKR) {
           touch[jj] = 1;
           R2=Reff*Reff;
@@ -391,6 +392,7 @@ void PairGranular::compute(int eflag, int vflag)
         } else {
           Fncrit = fabs(Fntot);
         }
+		Fscrit = tangential_coeffs[itype][jtype][2] * Fncrit;
 
         //------------------------------
         // tangential forces
@@ -446,7 +448,6 @@ void PairGranular::compute(int eflag, int vflag)
           fs3 = -k_tangential*history[2] - damp_tangential*vtr3;
 
           // rescale frictional displacements and forces if needed
-          Fscrit = tangential_coeffs[itype][jtype][2] * Fncrit;
           fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3);
           if (fs > Fscrit) {
             shrmag = sqrt(history[0]*history[0] + history[1]*history[1] +
@@ -464,8 +465,8 @@ void PairGranular::compute(int eflag, int vflag)
             } else fs1 = fs2 = fs3 = 0.0;
           }
         } else { // classic pair gran/hooke (no history)
-          fs = meff*damp_tangential*vrel;
+          fs = damp_tangential*vrel;
-          if (vrel != 0.0) Ft = MIN(Fne,fs) / vrel;
+          if (vrel != 0.0) Ft = MIN(Fscrit,fs) / vrel;
           else Ft = 0.0;
           fs1 = -Ft*vtr1;
           fs2 = -Ft*vtr2;
@@ -635,7 +636,7 @@ void PairGranular::compute(int eflag, int vflag)
             torque[j][2] -= torroll3;
           }
         }
-        if (evflag) ev_tally_xyz(i,j,nlocal,0,
+        if (evflag) ev_tally_xyz(i,j,nlocal,force->newton_pair,
             0.0,0.0,fx,fy,fz,delx,dely,delz);
       }
     }
@@ -1374,22 +1375,6 @@ double PairGranular::single(int i, int j, int itype, int jtype,
   vn2 = ny*vnnr;
   vn3 = nz*vnnr;
 
-  double *rmass = atom->rmass;
-  int *mask = atom->mask;
-  mi = rmass[i];
-  mj = rmass[j];
-  if (fix_rigid) {
-    if (mass_rigid[i] > 0.0) mi = mass_rigid[i];
-    if (mass_rigid[j] > 0.0) mj = mass_rigid[j];
-  }
-
-  meff = mi*mj / (mi+mj);
-  if (mask[i] & freeze_group_bit) meff = mj;
-  if (mask[j] & freeze_group_bit) meff = mi;
-
-  delta = radsum - r;
-  dR = delta*Reff;
-
   // tangential component
 
   vt1 = vr1 - vn1;
@@ -1407,6 +1392,9 @@ double PairGranular::single(int i, int j, int itype, int jtype,
   // if I or J part of rigid body, use body mass
   // if I or J is frozen, meff is other particle
 
+  double *rmass = atom->rmass;
+  int *mask = atom->mask;
+
   mi = rmass[i];
   mj = rmass[j];
   if (fix_rigid) {
@@ -1451,11 +1439,13 @@ double PairGranular::single(int i, int j, int itype, int jtype,
   }
 
   if (damping_model[itype][jtype] == VELOCITY) {
-    damp_normal = normal_coeffs[itype][jtype][1];
+    damp_normal = 1;
+  } else if (damping_model[itype][jtype] == MASS_VELOCITY) {
+    damp_normal = meff;
   } else if (damping_model[itype][jtype] == VISCOELASTIC) {
-    damp_normal = normal_coeffs[itype][jtype][1]*a*meff;
+    damp_normal = a*meff;
   } else if (damping_model[itype][jtype] == TSUJI) {
-    damp_normal = normal_coeffs[itype][jtype][1]*sqrt(meff*knfac);
+    damp_normal = sqrt(meff*knfac);
   }
 
   damp_normal_prefactor = normal_coeffs[itype][jtype][1]*damp_normal;
@@ -1506,6 +1496,7 @@ double PairGranular::single(int i, int j, int itype, int jtype,
   } else {
     Fncrit = fabs(Fntot);
   }
+  Fscrit = tangential_coeffs[itype][jtype][2] * Fncrit;
 
   //------------------------------
   // tangential forces
@@ -1518,13 +1509,6 @@ double PairGranular::single(int i, int j, int itype, int jtype,
       k_tangential *= a;
     } else if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) {
       k_tangential *= a;
-      // on unloading, rescale the shear displacements
-      if (a < history[3]) {
-        double factor = a/history[3];
-        history[0] *= factor;
-        history[1] *= factor;
-        history[2] *= factor;
-      }
     }
 
     shrmag = sqrt(history[0]*history[0] + history[1]*history[1] +
@@ -1535,35 +1519,34 @@ double PairGranular::single(int i, int j, int itype, int jtype,
     fs2 = -k_tangential*history[1] - damp_tangential*vtr2;
     fs3 = -k_tangential*history[2] - damp_tangential*vtr3;
 
-    // rescale frictional displacements and forces if needed
+    // rescale frictional forces if needed
-    Fscrit = tangential_coeffs[itype][jtype][2] * Fncrit;
     fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3);
     if (fs > Fscrit) {
       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);
         fs1 *= Fscrit/fs;
         fs2 *= Fscrit/fs;
         fs3 *= Fscrit/fs;
-      } else fs1 = fs2 = fs3 = 0.0;
+		fs *= Fscrit/fs;
+      } else fs1 = fs2 = fs3 = fs = 0.0;
     }
 
   // classic pair gran/hooke (no history)
   } else {
-    fs = meff*damp_tangential*vrel;
+    fs = damp_tangential*vrel;
-    if (vrel != 0.0) Ft = MIN(Fne,fs) / vrel;
+    if (vrel != 0.0) Ft = MIN(Fscrit,fs) / vrel;
     else Ft = 0.0;
     fs1 = -Ft*vtr1;
     fs2 = -Ft*vtr2;
     fs3 = -Ft*vtr3;
+	fs = Ft*vrel;
   }
 
   //****************************************
   // rolling resistance
   //****************************************
 
-  if (roll_model[itype][jtype] != ROLL_NONE) {
+  if ((roll_model[itype][jtype] != ROLL_NONE)
+      || (twist_model[itype][jtype] != TWIST_NONE)) {
     relrot1 = omega[i][0] - omega[j][0];
     relrot2 = omega[i][1] - omega[j][1];
     relrot3 = omega[i][2] - omega[j][2];
@@ -1601,7 +1584,8 @@ double PairGranular::single(int i, int j, int itype, int jtype,
         fr1 *= Frcrit/fr;
         fr2 *= Frcrit/fr;
         fr3 *= Frcrit/fr;
-      } else fr1 = fr2 = fr3 = 0.0;
+		fr *= Frcrit/fr;
+      } else fr1 = fr2 = fr3 = fr = 0.0;
     }
 
   }
@@ -1628,8 +1612,12 @@ double PairGranular::single(int i, int j, int itype, int jtype,
     Mtcrit = mu_twist*Fncrit; // critical torque (eq 44)
     if (fabs(magtortwist) > Mtcrit) {
       magtortwist = -Mtcrit * signtwist; // eq 34
+    } else magtortwist = 0.0;
   }
-  }
+	
+  // set force and return no energy
+	
+  fforce = Fntot*rinv;
 	
   // set single_extra quantities