Updates fix wall/gran to mirror recent updates to pair granular. Also some minor changes related to the limit_damping option

doc/src/pair_gran.rst

@@ -36,7 +36,7 @@ Syntax
 * xmu = static yield criterion (unitless value between 0.0 and 1.0e4)
 * dampflag = 0 or 1 if tangential damping force is excluded or included
 
-* keyword = *no_attraction*
+* keyword = *limit_damping*
 
   .. parsed-literal::
 
@@ -61,6 +61,8 @@ Examples
 
    pair_style gran/hooke/history 200000.0 NULL 50.0 NULL 0.5 1
    pair_style gran/hooke 200000.0 70000.0 50.0 30.0 0.5 0
+   pair_style gran/hooke 200000.0 70000.0 50.0 30.0 0.5 0 limit_damping
+
 
 Description
 """""""""""

doc/src/pair_granular.rst

@@ -24,7 +24,7 @@ Examples
    pair_coeff * * hooke 1000.0 50.0 tangential linear_history 500.0 1.0 0.4 damping mass_velocity
 
    pair_style granular
-   pair_coeff * * hertz 1000.0 50.0 tangential mindlin 1000.0 1.0 0.4
+   pair_coeff * * hertz 1000.0 50.0 tangential mindlin 1000.0 1.0 0.4 limit_damping
 
    pair_style granular
    pair_coeff * * hertz/material 1e8 0.3 0.3 tangential mindlin_rescale NULL 1.0 0.4 damping tsuji

src/GRANULAR/fix_wall_gran.cpp

@@ -53,7 +53,8 @@ enum{NONE,CONSTANT,EQUAL};
 enum {NORMAL_HOOKE, NORMAL_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};
 
@@ -216,7 +217,9 @@ FixWallGran::FixWallGran(LAMMPS *lmp, 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");
@@ -226,8 +229,14 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
             tangential_model = TANGENTIAL_MINDLIN;
           else if (strcmp(arg[iarg+1], "mindlin_rescale") == 0)
             tangential_model = TANGENTIAL_MINDLIN_RESCALE;
+          else if (strcmp(arg[iarg+1], "mindlin/force") == 0)
+            tangential_model = TANGENTIAL_MINDLIN_FORCE;
+          else if (strcmp(arg[iarg+1], "mindlin_rescale/force") == 0)
+            tangential_model = TANGENTIAL_MINDLIN_RESCALE_FORCE;
           if ((tangential_model == TANGENTIAL_MINDLIN ||
-               tangential_model == TANGENTIAL_MINDLIN_RESCALE) &&
+               tangential_model == TANGENTIAL_MINDLIN_RESCALE ||
+               tangential_model == TANGENTIAL_MINDLIN_FORCE ||
+               tangential_model == TANGENTIAL_MINDLIN_RESCALE_FORCE) &&
               (strcmp(arg[iarg+2], "NULL") == 0)) {
             if (normal_model == NORMAL_HERTZ || normal_model == NORMAL_HOOKE) {
               error->all(FLERR, "NULL setting for Mindlin tangential "
@@ -306,14 +315,20 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
       }
     }
     size_history = 3*tangential_history + 3*roll_history + twist_history;
+    //Unlike the pair style, the wall style does not have a 'touch'
+    //array. Hence, an additional entry in the history is used to
+    //determine if particles previously contacted for JKR cohesion purposes.
     if (normal_model == JKR) size_history += 1;
-    if (tangential_model == TANGENTIAL_MINDLIN_RESCALE) size_history += 1;
+    if (tangential_model == TANGENTIAL_MINDLIN_RESCALE ||
+        tangential_model == TANGENTIAL_MINDLIN_RESCALE_FORCE) size_history += 1;
   }
 
-  if(normal_model == JKR)
+  if (limit_damping and normal_model == JKR)
-    error->all(FLERR,"Cannot limit damping with JRK model");
+    error->all(FLERR,"Illegal pair_coeff command, "
-  if(normal_model == DMT)
+        "cannot limit damping with JRK model");
-    error->all(FLERR,"Cannot limit damping with DMT model");
+  if (limit_damping and normal_model == DMT)
+    error->all(FLERR,"Illegal pair_coeff command, "
+        "Cannot limit damping with DMT model");
 
   // wallstyle args
 
@@ -509,7 +524,8 @@ void FixWallGran::init()
     roll_history_index += 1;
     twist_history_index += 1;
   }
-  if (tangential_model == TANGENTIAL_MINDLIN_RESCALE) {
+  if (tangential_model == TANGENTIAL_MINDLIN_RESCALE ||
+      tangential_model == TANGENTIAL_MINDLIN_RESCALE_FORCE) {
     roll_history_index += 1;
     twist_history_index += 1;
   }
@@ -1120,7 +1136,8 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz,
   double signtwist, magtwist, magtortwist, Mtcrit;
   double tortwist1, tortwist2, tortwist3;
 
-  double shrmag,rsht;
+  double shrmag,rsht,prjmag;
+  bool frameupdate;
 
   r = sqrt(rsq);
   E = normal_coeffs[0];
@@ -1195,12 +1212,18 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz,
   Fdamp = -damp_normal_prefactor*vnnr;
 
   Fntot = Fne + Fdamp;
-  if(limit_damping and Fntot < 0.0) Fntot = 0.0;
+  if (limit_damping and Fntot < 0.0) Fntot = 0.0;
 
   //****************************************
   // 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;
@@ -1242,69 +1265,115 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz,
   int thist2 = thist1 + 1;
 
   if (tangential_history) {
-    if (tangential_model == TANGENTIAL_MINDLIN) {
+    if (tangential_model == TANGENTIAL_MINDLIN ||
+        tangential_model == TANGENTIAL_MINDLIN_FORCE) {
       k_tangential *= a;
     }
-    else if (tangential_model == TANGENTIAL_MINDLIN_RESCALE) {
+    else if (tangential_model ==
+             TANGENTIAL_MINDLIN_RESCALE ||
+             tangential_model ==
+             TANGENTIAL_MINDLIN_RESCALE_FORCE){
       k_tangential *= a;
-      if (a < history[3]) { //On unloading, rescale the shear displacements
+      // on unloading, rescale the shear displacements/force
+      if (a < history[thist2+1]) {
         double factor = a/history[thist2+1];
         history[thist0] *= factor;
         history[thist1] *= factor;
         history[thist2] *= factor;
       }
     }
-    shrmag = sqrt(history[thist0]*history[thist0] +
+
-                  history[thist1]*history[thist1] +
-                  history[thist2]*history[thist2]);
 
     // rotate and update displacements.
     // see e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235
     if (history_update) {
       rsht = history[thist0]*nx + history[thist1]*ny + history[thist2]*nz;
-      if (fabs(rsht) < EPSILON) rsht = 0;
+      if (tangential_model == TANGENTIAL_MINDLIN_FORCE ||
-      if (rsht > 0) {
+          tangential_model == TANGENTIAL_MINDLIN_RESCALE_FORCE)
-        // if rhst == shrmag, contacting pair has rotated 90 deg in one step,
+        frameupdate = fabs(rsht) > EPSILON*Fscrit;
-        // in which case you deserve a crash!
+      else
-        scalefac = shrmag/(shrmag - rsht);
+        frameupdate = fabs(rsht)*k_tangential > EPSILON*Fscrit;
+      if (frameupdate) {
+        shrmag = sqrt(history[thist0]*history[thist0] +
+            history[thist1]*history[thist1] +
+            history[thist2]*history[thist2]);
+        // projection
         history[thist0] -= rsht*nx;
         history[thist1] -= rsht*ny;
         history[thist2] -= 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[thist0] *= scalefac;
         history[thist1] *= scalefac;
         history[thist2] *= scalefac;
       }
       // update history
+      if (tangential_model == TANGENTIAL_HISTORY ||
+          tangential_model == TANGENTIAL_MINDLIN ||
+          tangential_model == TANGENTIAL_MINDLIN_RESCALE) {
       history[thist0] += vtr1*dt;
       history[thist1] += vtr2*dt;
       history[thist2] += vtr3*dt;
+      } else{
+        // tangential force
+        // see e.g. eq. 18 of Thornton et al, Pow. Tech. 2013, v223,p30-46
+        history[thist0] -= k_tangential*vtr1*dt;
+        history[thist1] -= k_tangential*vtr2*dt;
+        history[thist2] -= k_tangential*vtr3*dt;
+      }
+      if (tangential_model == TANGENTIAL_MINDLIN_RESCALE ||
+          tangential_model == TANGENTIAL_MINDLIN_RESCALE_FORCE)
+        history[thist2+1] = a;
     }
 
     // tangential forces = history + tangential velocity damping
-    fs1 = -k_tangential*history[thist0] - damp_tangential*vtr1;
+    if (tangential_model == TANGENTIAL_HISTORY ||
-    fs2 = -k_tangential*history[thist1] - damp_tangential*vtr2;
+        tangential_model == TANGENTIAL_MINDLIN ||
-    fs3 = -k_tangential*history[thist2] - damp_tangential*vtr3;
+        tangential_model == TANGENTIAL_MINDLIN_RESCALE) {
+      fs1 = -k_tangential*history[thist0] - damp_tangential*vtr1;
+      fs2 = -k_tangential*history[thist1] - damp_tangential*vtr2;
+      fs3 = -k_tangential*history[thist2] - damp_tangential*vtr3;
+    } else {
+      fs1 = history[thist0] - damp_tangential*vtr1;
+      fs2 = history[thist1] - damp_tangential*vtr2;
+      fs3 = history[thist2] - damp_tangential*vtr3;
+    }
 
     // rescale frictional displacements and forces if needed
     Fscrit = tangential_coeffs[2] * Fncrit;
     fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3);
     if (fs > Fscrit) {
+      shrmag = sqrt(history[thist0]*history[thist0] +
+          history[thist1]*history[thist1] +
+          history[thist2]*history[thist2]);
       if (shrmag != 0.0) {
-        history[thist0] = -1.0/k_tangential*(Fscrit*fs1/fs +
+        if (tangential_model == TANGENTIAL_HISTORY ||
-                                             damp_tangential*vtr1);
+            tangential_model == TANGENTIAL_MINDLIN ||
-        history[thist1] = -1.0/k_tangential*(Fscrit*fs2/fs +
+            tangential_model ==
-                                             damp_tangential*vtr2);
+                TANGENTIAL_MINDLIN_RESCALE) {
-        history[thist2] = -1.0/k_tangential*(Fscrit*fs3/fs +
+          history[thist0] = -1.0/k_tangential*(Fscrit*fs1/fs +
-                                             damp_tangential*vtr3);
+              damp_tangential*vtr1);
+          history[thist1] = -1.0/k_tangential*(Fscrit*fs2/fs +
+              damp_tangential*vtr2);
+          history[thist2] = -1.0/k_tangential*(Fscrit*fs3/fs +
+              damp_tangential*vtr3);
+        } else {
+          history[thist0] = Fscrit*fs1/fs + damp_tangential*vtr1;
+          history[thist1] = Fscrit*fs2/fs + damp_tangential*vtr2;
+          history[thist2] = Fscrit*fs3/fs + damp_tangential*vtr3;
+        }
         fs1 *= Fscrit/fs;
         fs2 *= Fscrit/fs;
         fs3 *= Fscrit/fs;
       } 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;
@@ -1315,14 +1384,14 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz,
   // rolling resistance
   //****************************************
 
-  if (roll_model != ROLL_NONE || twist_model != NONE) {
+  if (roll_model != ROLL_NONE || twist_model != TWIST_NONE) {
     relrot1 = omega[0];
     relrot2 = omega[1];
     relrot3 = omega[2];
   }
   if (roll_model != ROLL_NONE) {
-
+    // rolling velocity,
-    // rolling velocity, see eq. 31 of Wang et al, Particuology v 23, p 49 (2015)
+    // see eq. 31 of Wang et al, Particuology v 23, p 49 (2015)
     // This is different from the Marshall papers,
     // which use the Bagi/Kuhn formulation
     // for rolling velocity (see Wang et al for why the latter is wrong)
@@ -1334,21 +1403,29 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz,
     int rhist1 = rhist0 + 1;
     int rhist2 = rhist1 + 1;
 
-    // rolling displacement
+    k_roll = roll_coeffs[0];
-    rollmag = sqrt(history[rhist0]*history[rhist0] +
+    damp_roll = roll_coeffs[1];
-                   history[rhist1]*history[rhist1] +
+    Frcrit = roll_coeffs[2] * Fncrit;
-                   history[rhist2]*history[rhist2]);
-
-    rolldotn = history[rhist0]*nx + history[rhist1]*ny + history[rhist2]*nz;
 
     if (history_update) {
-      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;
-        scalefac = rollmag/(rollmag - rolldotn);
+      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;
@@ -1358,17 +1435,16 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz,
       history[rhist2] += vrl3*dt;
     }
 
-    k_roll = roll_coeffs[0];
-    damp_roll = roll_coeffs[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[2] * Fncrit;
-
     fr = sqrt(fr1*fr1 + fr2*fr2 + fr3*fr3);
     if (fr > Frcrit) {
+      rollmag = sqrt(history[rhist0]*history[rhist0] +
+          history[rhist1]*history[rhist1] +
+          history[rhist2]*history[rhist2]);
       if (rollmag != 0.0) {
         history[rhist0] = -1.0/k_roll*(Frcrit*fr1/fr + damp_roll*vrl1);
         history[rhist1] = -1.0/k_roll*(Frcrit*fr2/fr + damp_roll*vrl2);

src/GRANULAR/pair_gran_hooke_history.cpp

@@ -237,7 +237,7 @@ void PairGranHookeHistory::compute(int eflag, int vflag)
 
         damp = meff*gamman*vnnr*rsqinv;
         ccel = kn*(radsum-r)*rinv - damp;
-        if(limit_damping and ccel < 0.0) ccel = 0.0;
+        if (limit_damping and ccel < 0.0) ccel = 0.0;
 
         // relative velocities
 
@@ -372,8 +372,8 @@ void PairGranHookeHistory::settings(int narg, char **arg)
   if (dampflag == 0) gammat = 0.0;
   
   limit_damping = 0;
-  if(narg == 7) {
+  if (narg == 7) {
-    if(strcmp(arg[6], "limit_damping") == 0) limit_damping = 1;
+    if (strcmp(arg[6], "limit_damping") == 0) limit_damping = 1;
     else error->all(FLERR,"Illegal pair_style command");
   }
 

src/GRANULAR/pair_granular.cpp

@@ -365,13 +365,13 @@ void PairGranular::compute(int eflag, int vflag)
           damp_normal = a*meff;
         } else if (damping_model[itype][jtype] == TSUJI) {
           damp_normal = sqrt(meff*knfac);
-        }
+        } else damp_normal = 0.0;
 
         damp_normal_prefactor = normal_coeffs[itype][jtype][1]*damp_normal;
         Fdamp = -damp_normal_prefactor*vnnr;
 
         Fntot = Fne + Fdamp;
-        if(limit_damping[itype][jtype] and Fntot < 0.0) Fntot = 0.0;
+        if (limit_damping[itype][jtype] and Fntot < 0.0) Fntot = 0.0;
 
         //****************************************
         // tangential force, including history effects
@@ -542,20 +542,17 @@ void PairGranular::compute(int eflag, int vflag)
           relrot1 = omega[i][0] - omega[j][0];
           relrot2 = omega[i][1] - omega[j][1];
           relrot3 = omega[i][2] - omega[j][2];
-          // rolling velocity,
-          // see eq. 31 of Wang et al, Particuology v 23, p 49 (2015)
-          // this is different from the Marshall papers,
-          // which use the Bagi/Kuhn formulation
-          // for rolling velocity (see Wang et al for why the latter is wrong)
-          // - 0.5*((radj-radi)/radsum)*vtr1;
-          // - 0.5*((radj-radi)/radsum)*vtr2;
-          // - 0.5*((radj-radi)/radsum)*vtr3;
         }
         //****************************************
         // rolling resistance
         //****************************************
 
         if (roll_model[itype][jtype] != ROLL_NONE) {
+          // rolling velocity,
+          // see eq. 31 of Wang et al, Particuology v 23, p 49 (2015)
+          // this is different from the Marshall papers,
+          // which use the Bagi/Kuhn formulation
+          // for rolling velocity (see Wang et al for why the latter is wrong)
           vrl1 = Reff*(relrot2*nz - relrot3*ny);
           vrl2 = Reff*(relrot3*nx - relrot1*nz);
           vrl3 = Reff*(relrot1*ny - relrot2*nx);
@@ -974,12 +971,12 @@ void PairGranular::coeff(int narg, char **arg)
     } else if (strcmp(arg[iarg], "limit_damping") == 0) {
       ld_flag = 1;        
       iarg += 1;
-    } else error->all(FLERR, "Illegal pair coeff command");
+    } else error->all(FLERR, "Illegal pair_coeff command");
   }
 
   // error not to specify normal or tangential model
   if ((normal_model_one < 0) || (tangential_model_one < 0))
-    error->all(FLERR, "Illegal pair coeff command, "
+    error->all(FLERR, "Illegal pair_coeff command, "
                "must specify normal or tangential contact model");
 
   int count = 0;
@@ -991,11 +988,13 @@ void PairGranular::coeff(int narg, char **arg)
         27.467*powint(cor,4)-18.022*powint(cor,5)+4.8218*powint(cor,6);
   } else damp = normal_coeffs_one[1];
 
-  if(ld_flag and normal_model_one == JKR)
+  if (ld_flag and normal_model_one == JKR)
-    error->all(FLERR,"Cannot limit damping with JKR model");
+    error->all(FLERR,"Illegal pair_coeff command, "
+        "Cannot limit damping with JKR model");
 
-  if(ld_flag and normal_model_one == DMT)
+  if (ld_flag and normal_model_one == DMT)
-    error->all(FLERR,"Cannot limit damping with DMT model");
+    error->all(FLERR,"Illegal pair_coeff command, "
+        "Cannot limit damping with DMT model");
 
   for (int i = ilo; i <= ihi; i++) {
     for (int j = MAX(jlo,i); j <= jhi; j++) {