Commit 2cc35ee7 authored by Dan S. Bolintineanu's avatar Dan S. Bolintineanu
Browse files

A few additional enhancements to pair granular and fix wall granular option: 

- NULL option for tangential stiffness, to set it based on shear modulus
- calculation of effective shear moduli from elastic moduli and Poisson's ratio
- updates to doc page example syntax
parent e0c935b5

doc/src/fix_wall_gran.txt

+2 −2
Original line number Diff line number Diff line
@@ -49,8 +49,8 @@ fix 1 all wall/gran hooke 200000.0 NULL 50.0 NULL 0.5 0 xplane -10.0 10.0 
fix 1 all wall/gran hooke/history 200000.0 NULL 50.0 NULL 0.5 0 zplane 0.0 NULL

fix 2 all wall/gran hooke 100000.0 20000.0 50.0 30.0 0.5 1 zcylinder 15.0 wiggle z 3.0 2.0 

fix 3 all wall/gran granular hooke 1000.0 50.0 tangential linear_nohistory 1.0 0.4 zplane 0.0 NULL

fix 4 all wall/gran granular jkr 1000.0 50.0 tangential linear_history 800.0 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall zcylinder 15.0 wiggle z 3.0 2.0

fix 5 all wall/gran granular dmt 1000.0 50.0 0.3 10.0 tangential linear_history 800.0 0.5 0.1 roll sds 500.0 200.0 0.1 twisting marshall zplane 0.0 NULL :pre

fix 4 all wall/gran granular jkr 1000.0 50.0 0.3 5.0 tangential mindlin 800.0 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall zcylinder 15.0 wiggle z 3.0 2.0

fix 5 all wall/gran granular dmt 1000.0 50.0 0.3 10.0 tangential mindlin 800.0 0.5 0.1 roll sds 500.0 200.0 0.1 twisting marshall zplane 0.0 NULL :pre



[Description:]

doc/src/pair_granular.txt

+24 −15
Original line number Diff line number Diff line
@@ -28,10 +28,10 @@ pair_style granular 
pair_coeff * * hooke 1000.0 50.0 tangential linear_nohistory 1.0 0.4 :pre



pair_style granular

pair_coeff * * hertz 1000.0 50.0 tangential mindlin 800.0 1.0 0.4 :pre

pair_coeff * * hertz 1000.0 50.0 tangential mindlin NULL 1.0 0.4 :pre



pair_style granular

pair_coeff * * hertz 1000.0 0.3 tangential linear_history 800.0 1.0 0.4 damping tsuji :pre

pair_coeff * * hertz/material 1e8 0.3 tangential mindlin_rescale NULL 1.0 0.4 damping tsuji :pre



pair_style granular

pair_coeff 1 1 jkr 1000.0 50.0 tangential mindlin 800.0 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall
@@ -39,8 +39,8 @@ pair_coeff 2 2 hertz 200.0 20.0 tangential linear_history 300.0 1.0 0.1 rolling 


pair_style granular

pair_coeff 1 1 hertz 1000.0 50.0 tangential mindlin 800.0 0.5 0.5 rolling sds 500.0 200.0 0.5 twisting marshall 

pair_coeff 2 2 dmt 1000.0 50.0 0.3 10.0 tangential linear_history 800.0 0.5 0.1 roll sds 500.0 200.0 0.1 twisting marshall

pair_coeff 1 2 dmt 1000.0 50.0 0.3 10.0 tangential linear_history 800.0 0.5 0.1 roll sds 500.0 200.0 0.1 twisting marshall :pre

pair_coeff 2 2 dmt 1000.0 50.0 0.3 10.0 tangential mindlin 800.0 0.5 0.1 roll sds 500.0 200.0 0.1 twisting marshall

pair_coeff 1 2 dmt 1000.0 50.0 0.3 10.0 tangential mindlin 800.0 0.5 0.1 roll sds 500.0 200.0 0.1 twisting marshall :pre



[Description:]
@@ -220,8 +220,8 @@ expected parameters are as follows: 
 

{linear_nohistory} : \(x_\{\gamma,t\}\), \(\mu_s\)

{linear_history} : \(k_t\), \(x_\{\gamma,t\}\), \(\mu_s\) 

{mindlin} : \(k_t\), \(x_\{\gamma,t\}\), \(\mu_s\) 

{mindlin_rescale} : \(k_t\), \(x_\{\gamma,t\}\), \(\mu_s\) :ol

{mindlin} : \(k_t\) or NULL, \(x_\{\gamma,t\}\), \(\mu_s\) 

{mindlin_rescale} : \(k_t\) or NULL, \(x_\{\gamma,t\}\), \(\mu_s\) :ol



Here, \(x_\{\gamma,t\}\) is a dimensionless multiplier for the normal damping \(\eta_n\)

that determines the magnitude of the
@@ -337,7 +337,7 @@ to induce a torque on each particle according to: 
\mathbf\{\tau\}_j = -(R_j - 0.5 \delta) \mathbf\{n\} \times \mathbf\{F\}_t

\end\{equation\}



For {tangential mindlin}, the Mindlin no-slip solution is used, which differs from the {linear_history}

For {tangential mindlin}, the "Mindlin"_#Mindlin1949 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:



\begin\{equation\}
@@ -346,18 +346,27 @@ option by an additional factor of {a}, the radius of the contact region. The tan 


Here, {a} is the radius of the contact region, given by \(a = \delta R\) 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.

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 \(E\) by \(G = E/(2(1+\nu))\), where \(\nu\) 

is Poisson's ratio. This can also be achieved by specifying {NULL} for \(k_t\), in which case

a normal contact model that specifies material parameters \(E\) and \(\nu\) is required (e.g. {hertz/material},

{dmt} or {jkr}). In this case, mixing of shear moduli for different particle types {i} and {j} is done according

to:

\begin\{equation\}

1/G = 2(2-\nu_i)(1+\nu_i)/E_i + 2(2-\nu_j)(1+\nu_j)/E_j 

\end\{equation\}



The {mindlin_rescale} option uses the same form as {mindlin}, but the magnitude of the tangential 

displacement is re-scaled as the contact unloads, i.e. if \(a < a_\{t_n-1\}\):

displacement is re-scaled as the contact unloads, i.e. if \(a < a_\{t_\{n-1\}\}\):

\begin\{equation\}

\mathbf\{\xi\} = \mathbf\{xi_\{t_n-1\}\} \frac\{a\}\{a_\{t_n-1\}\}

\mathbf\{\xi\} = \mathbf\{\xi_\{t_\{n-1\}\}\} \frac\{a\}\{a_\{t_\{n-1\}\}\}

\end\{equation\}



This accounts for the fact that a decrease in the contact area upon unloading leads to the contact

Here, \(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 ("Walton"_#WaltonPC). See also discussion in "Thornton et al, 2013"_#Thornton2013,

particularly equation 18(b) of that work and associated discussion.

without the rescaling above ("Walton"_#WaltonPC ). See also discussion in "Thornton et al, 2013"_#Thornton2013

, particularly equation 18(b) of that work and associated discussion.

 

:line
@@ -649,10 +658,10 @@ J. Phys. D: Appl. Phys. 24 1942 
[(Mindlin, 1949)] Mindlin, R. D. (1949). Compliance of elastic bodies in contact.

J. Appl. Mech., ASME 16, 259-268.



:line(Thornton2013)

:link(Thornton2013)

[(Thornton et al, 2013)] Thornton, C., Cummins, S. J., & Cleary, P. W. (2013). 

An investigation of the comparative behaviour of alternative contact force models 

during inelastic collisions. Powder Technology, 233, 30-46.



:line(WaltonPC)

:link(WaltonPC)

[(Otis R. Walton)] Walton, O.R., Personal Communication

src/GRANULAR/fix_wall_gran.cpp

+38 −5
Original line number Diff line number Diff line
@@ -54,7 +54,7 @@ enum{NONE,CONSTANT,EQUAL}; 


enum {NORMAL_HOOKE, NORMAL_HERTZ, HERTZ_MATERIAL, DMT, JKR};

enum {VELOCITY, VISCOELASTIC, TSUJI};

enum {TANGENTIAL_NOHISTORY, TANGENTIAL_HISTORY, TANGENTIAL_MINDLIN};

enum {TANGENTIAL_NOHISTORY, TANGENTIAL_HISTORY, TANGENTIAL_MINDLIN, TANGENTIAL_MINDLIN_RESCALE};

enum {TWIST_NONE, TWIST_SDS, TWIST_MARSHALL};

enum {ROLL_NONE, ROLL_SDS};
@@ -198,11 +198,24 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : 
          tangential_coeffs[2] = force->numeric(FLERR,arg[iarg+3]); //friction coeff.

          iarg += 4;

        }

        else if (strcmp(arg[iarg+1], "linear_history") == 0){

        else if ((strcmp(arg[iarg+1], "linear_history") == 0) ||

            (strcmp(arg[iarg+1], "mindlin") == 0) ||

            (strcmp(arg[iarg+1], "mindlin_rescale") == 0)){

          if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for tangential model");

          tangential_model = TANGENTIAL_HISTORY;

          tangential_history = 1;

          if (strcmp(arg[iarg+1], "linear_history") == 0) tangential_model = TANGENTIAL_HISTORY;

          else if (strcmp(arg[iarg+1], "mindlin") == 0) tangential_model = TANGENTIAL_MINDLIN;

          else if (strcmp(arg[iarg+1], "mindlin_rescale") == 0) tangential_model = TANGENTIAL_MINDLIN_RESCALE;

          if ((tangential_model == TANGENTIAL_MINDLIN || tangential_model == TANGENTIAL_MINDLIN_RESCALE) &&

              (strcmp(arg[iarg+2], "NULL") == 0)){

            if (normal_model == NORMAL_HERTZ || normal_model == NORMAL_HOOKE){

              error->all(FLERR, "NULL setting for Mindlin tangential stiffness requires a normal contact model that specifies material properties");

            }

            tangential_coeffs[0] = 4*(2-poiss)*(1+poiss)/Emod;

          }

          else{

            tangential_coeffs[0] = force->numeric(FLERR,arg[iarg+2]); //kt

          }

          tangential_history = 1;

          tangential_coeffs[1] = force->numeric(FLERR,arg[iarg+3]); //gammat

          tangential_coeffs[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff.

          iarg += 5;
@@ -265,7 +278,9 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : 
        error->all(FLERR, "Illegal fix wall/gran command");

      }

    }

    size_history = (normal_model == JKR) + 3*tangential_history + 3*roll_history + twist_history;

    size_history = 3*tangential_history + 3*roll_history + twist_history;

    if (normal_model == JKR) size_history += 1;

    if (tangential_model == TANGENTIAL_MINDLIN_RESCALE) size_history += 1;

  }



  // wallstyle args
@@ -466,6 +481,10 @@ void FixWallGran::init() 
    roll_history_index += 1;

    twist_history_index += 1;

  }

  if (tangential_model == TANGENTIAL_MINDLIN_RESCALE){

    roll_history_index += 1;

    twist_history_index += 1;

  }



  if (damping_model == TSUJI){

    double cor = normal_coeffs[1];
@@ -473,6 +492,8 @@ void FixWallGran::init() 
        27.467*pow(cor,4)-18.022*pow(cor,5)+

        4.8218*pow(cor,6);

  }





}



/* ---------------------------------------------------------------------- */
@@ -1177,6 +1198,18 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz, 
 int thist2 = thist1 + 1;



 if (tangential_history){

   if (tangential_model == TANGENTIAL_MINDLIN){

     k_tangential *= a;

   }

   else if (tangential_model == TANGENTIAL_MINDLIN_RESCALE){

     k_tangential *= a;

     if (a < history[3]){ //On unloading, rescale the shear displacements

       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]);

src/GRANULAR/pair_granular.cpp

+47 −13
Original line number Diff line number Diff line
@@ -211,8 +211,10 @@ void PairGranular::compute(int eflag, int vflag) 
  ilist = list->ilist;

  numneigh = list->numneigh;

  firstneigh = list->firstneigh;

  if (use_history){

    firsttouch = fix_history->firstflag;

    firsthistory = fix_history->firstvalue;

  }



  for (ii = 0; ii < inum; ii++) {

    i = ilist[ii];
@@ -222,8 +224,10 @@ void PairGranular::compute(int eflag, int vflag) 
    ztmp = x[i][2];

    itype = type[i];

    radi = radius[i];

    if (use_history){

      touch = firsttouch[i];

      allhistory = firsthistory[i];

    }

    jlist = firstneigh[i];

    jnum = numneigh[i];
@@ -263,10 +267,12 @@ void PairGranular::compute(int eflag, int vflag) 


      if (!touchflag){

        // unset non-touching neighbors

        if (use_history){

          touch[jj] = 0;

          history = &allhistory[size_history*jj];

          for (int k = 0; k < size_history; k++) history[k] = 0.0;

        }

      }

      else{

        r = sqrt(rsq);

        rinv = 1.0/r;
@@ -788,7 +794,16 @@ void PairGranular::coeff(int narg, char **arg) 
        else if (strcmp(arg[iarg+1], "mindlin") == 0) tangential_model_one = TANGENTIAL_MINDLIN;

        else if (strcmp(arg[iarg+1], "mindlin_rescale") == 0) tangential_model_one = TANGENTIAL_MINDLIN_RESCALE;

        tangential_history = 1;

        if ((tangential_model_one == TANGENTIAL_MINDLIN || tangential_model_one == TANGENTIAL_MINDLIN_RESCALE) &&

            (strcmp(arg[iarg+2], "NULL") == 0)){

          if (normal_model_one == HERTZ || normal_model_one == HOOKE){

            error->all(FLERR, "NULL setting for Mindlin tangential stiffness requires a normal contact model that specifies material properties");

          }

          tangential_coeffs_one[0] = -1;

        }

        else{

          tangential_coeffs_one[0] = force->numeric(FLERR,arg[iarg+2]); //kt

        }

        tangential_coeffs_one[1] = force->numeric(FLERR,arg[iarg+3]); //gammat

        tangential_coeffs_one[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff.

        iarg += 5;
@@ -879,7 +894,13 @@ void PairGranular::coeff(int narg, char **arg) 
      damping_model[i][j] = damping_model[j][i] = damping_model_one;



      tangential_model[i][j] = tangential_model[j][i] = tangential_model_one;

      for (int k = 0; k < 3; k++)

      if (tangential_coeffs_one[0] == -1){

        tangential_coeffs[i][j][0] = tangential_coeffs[j][i][0] = 8*mix_stiffnessG(Emod[i][j], Emod[i][j], poiss[i][j], poiss[i][j]);

      }

      else{

        tangential_coeffs[i][j][0] = tangential_coeffs[j][i][0] = tangential_coeffs_one[0];

      }

      for (int k = 1; k < 3; k++)

        tangential_coeffs[i][j][k] = tangential_coeffs[j][i][k] = tangential_coeffs_one[k];



      roll_model[i][j] = roll_model[j][i] = roll_model_one;
@@ -1276,7 +1297,7 @@ double PairGranular::single(int i, int j, int itype, int jtype, 
    touchflag = (rsq <= radsum*radsum);

  }



  if (touchflag){

  if (!touchflag){

    fforce = 0.0;

    for (int m = 0; m < single_extra; m++) svector[m] = 0.0;

    return 0.0;
@@ -1373,8 +1394,8 @@ double PairGranular::single(int i, int j, int itype, int jtype, 
  }

  else{

    knfac = E;

    a = sqrt(dR);

    Fne = knfac*delta;

    a = sqrt(dR);

    if (normal_model[itype][jtype] != HOOKE){

      Fne *= a;

      knfac *= a;
@@ -1451,6 +1472,19 @@ 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){

      k_tangential *= a;

      if (a < history[3]){ //On unloading, rescale the shear displacements

        double factor = a/history[3];

        history[0] *= factor;

        history[1] *= factor;

        history[2] *= factor;

      }

    }



    shrmag = sqrt(history[0]*history[0] + history[1]*history[1] +

        history[2]*history[2]);
@@ -1617,9 +1651,9 @@ double PairGranular::mix_stiffnessE(double Eii, double Ejj, double poisii, doubl 
	 mixing of shear modulus (G)

------------------------------------------------------------------------ */



double PairGranular::mix_stiffnessG(double Gii, double Gjj, double poisii, double poisjj)

double PairGranular::mix_stiffnessG(double Eii, double Ejj, double poisii, double poisjj)

{

  return 1/((2.0 -poisii)/Gii+(2.0-poisjj)/Gjj);

  return 1/((2*(2-poisii)*(1+poisii)/Eii) + (2*(2-poisjj)*(1+poisjj)/Ejj));

}



/* ----------------------------------------------------------------------

CRA Git | Maintained and supported by SUSTech CRA and CCSE