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

Additional changes to pair_granular: 

-added mindlin and mindlin/rescale tangential model options
-torque from tangential force is now applied at the same point on both contacting particles
-updated doc page to reflect changes above
parent 4ee98d18

doc/src/pair_granular.txt

+48 −8
Original line number Diff line number Diff line
@@ -28,17 +28,17 @@ 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 linear_history 800.0 1.0 0.4 :pre

pair_coeff * * hertz 1000.0 50.0 tangential mindlin 800.0 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_style granular

pair_coeff 1 1 jkr 1000.0 50.0 tangential linear_history 800.0 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall 

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 

pair_coeff 2 2 hertz 200.0 20.0 tangential linear_history 300.0 1.0 0.1 rolling sds 200.0 100.0 0.1 twisting marshall :pre



pair_style granular

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

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
@@ -219,7 +219,9 @@ choice and associated parameters. Currently supported tangential model choices a 
expected parameters are as follows:

 

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

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

{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



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

that determines the magnitude of the
@@ -273,6 +275,10 @@ F_\{n0\} = \|\mathbf\{F\}_ne + 2 F_\{pulloff\}\| 


Where \(F_\{pulloff\} = 3\pi \gamma R \) for {jkr}, and \(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.



For {tangential linear_history}, the tangential force is given by:



\begin\{equation\}
@@ -320,17 +326,39 @@ equation 20 and related discussion): 
\end\{equation\}



The tangential force is added to the total normal force (elastic plus damping) to produce the total force

on the particle. The tangential force also acts at the contact point to induce a torque on each 

particle according to:

on the particle. The tangential force also acts at the contact point (defined as the center of the overlap region)

to induce a torque on each particle according to:



\begin\{equation\}

\mathbf\{\tau\}_i = -R_i \mathbf\{n\} \times \mathbf\{F\}_t

\mathbf\{\tau\}_i = -(R_i - 0.5 \delta) \mathbf\{n\} \times \mathbf\{F\}_t

\end\{equation\}



\begin\{equation\}

\mathbf\{\tau\}_j = -R_j \mathbf\{n\} \times \mathbf\{F\}_t

\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}

option by an additional factor of {a}, the radius of the contact region. The tangential force is given by:



\begin\{equation\}                                                                                                              

\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 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.



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\}\):

\begin\{equation\}

\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

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.

 

:line



The optional {rolling} keyword enables rolling friction, which resists pure rolling
@@ -616,3 +644,15 @@ Rolling and sliding in 3-D discrete element models. Particuology, 23, 49-55. 
:link(Thornton1991)

[(Thornton, 1991)] Thornton, C. (1991). Interparticle sliding in the presence of adhesion.

J. Phys. D: Appl. Phys. 24 1942



:link(Mindlin1949)

[(Mindlin, 1949)] Mindlin, R. D. (1949). Compliance of elastic bodies in contact.

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



:line(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)

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

src/GRANULAR/pair_granular.cpp

+66 −22
Original line number Diff line number Diff line
@@ -52,7 +52,7 @@ using namespace MathConst; 


enum {HOOKE, 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};
@@ -77,6 +77,8 @@ PairGranular::PairGranular(LAMMPS *lmp) : Pair(lmp) 
  maxrad_dynamic = NULL;

  maxrad_frozen = NULL;



  history_transfer_factors = NULL;



  dt = update->dt;



  // set comm size needed by this Pair if used with fix rigid
@@ -131,7 +133,7 @@ void PairGranular::compute(int eflag, int vflag) 
  int i,j,ii,jj,inum,jnum,itype,jtype;

  double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz,nx,ny,nz;

  double radi,radj,radsum,rsq,r,rinv;

  double Reff, delta, dR, dR2;

  double Reff, delta, dR, dR2, dist_to_contact;



  double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;

  double wr1,wr2,wr3;
@@ -397,16 +399,27 @@ void PairGranular::compute(int eflag, int vflag) 
        damp_tangential = tangential_coeffs[itype][jtype][1]*damp_normal_prefactor;



        if (tangential_history){

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

              history[2]*history[2]);



          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;

            }

          }

          // Rotate and update displacements.

          // 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){

              scalefac = shrmag/(shrmag - rsht); //if rhst == shrmag, contacting pair has rotated 90 deg. in one step, in which case you deserve a crash!

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

                                               history[2]*history[2]);

              scalefac = shrmag/(shrmag - rsht); //if rsht == shrmag, contacting pair has rotated 90 deg. in one step, in which case you deserve a crash!

              history[0] -= rsht*nx;

              history[1] -= rsht*ny;

              history[2] -= rsht*nz;
@@ -419,6 +432,7 @@ void PairGranular::compute(int eflag, int vflag) 
            history[0] += vtr1*dt;

            history[1] += vtr2*dt;

            history[2] += vtr3*dt;

            if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) history[3] = a;

          }



          // tangential forces = history + tangential velocity damping
@@ -430,6 +444,8 @@ void PairGranular::compute(int eflag, int vflag) 
          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] +

                                    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);
@@ -469,16 +485,13 @@ void PairGranular::compute(int eflag, int vflag) 
          int rhist1 = rhist0 + 1;

          int rhist2 = rhist1 + 1;



          // Rolling displacement

          rollmag = sqrt(history[rhist0]*history[rhist0] +

              history[rhist1]*history[rhist1] +

              history[rhist2]*history[rhist2]);



          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;
@@ -504,6 +517,9 @@ void PairGranular::compute(int eflag, int vflag) 


          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);
@@ -555,9 +571,10 @@ void PairGranular::compute(int eflag, int vflag) 
        tor2 = nz*fs1 - nx*fs3;

        tor3 = nx*fs2 - ny*fs1;



        torque[i][0] -= radi*tor1;

        torque[i][1] -= radi*tor2;

        torque[i][2] -= radi*tor3;

	dist_to_contact = radi-0.5*delta;

        torque[i][0] -= dist_to_contact*tor1;

        torque[i][1] -= dist_to_contact*tor2;

        torque[i][2] -= dist_to_contact*tor3;



        if (twist_model[itype][jtype] != TWIST_NONE){

          tortwist1 = magtortwist * nx;
@@ -584,9 +601,10 @@ void PairGranular::compute(int eflag, int vflag) 
          f[j][1] -= fy;

          f[j][2] -= fz;



          torque[j][0] -= radj*tor1;

          torque[j][1] -= radj*tor2;

          torque[j][2] -= radj*tor3;

	  dist_to_contact = radj-0.5*delta;

          torque[j][0] -= dist_to_contact*tor1;

          torque[j][1] -= dist_to_contact*tor2;

          torque[j][2] -= dist_to_contact*tor3;



          if (twist_model[itype][jtype] != TWIST_NONE){

            torque[j][0] -= tortwist1;
@@ -762,9 +780,13 @@ void PairGranular::coeff(int narg, char **arg) 
        tangential_coeffs_one[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_one = TANGENTIAL_HISTORY;

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

        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;

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

        tangential_coeffs_one[1] = force->numeric(FLERR,arg[iarg+3]); //gammat
@@ -896,11 +918,11 @@ void PairGranular::init_style() 
    error->all(FLERR,"Pair granular requires ghost atoms store velocity");



  // Determine whether we need a granular neigh list, how large it needs to be

  use_history = tangential_history || roll_history || twist_history;

  use_history = normal_history || tangential_history || roll_history || twist_history;



  //For JKR, will need fix/neigh/history to keep track of touch arrays

  for (int i = 1; i <= atom->ntypes; i++)

    for (int j = 1; j <= atom->ntypes; j++)

    for (int j = i; j <= atom->ntypes; j++)

      if (normal_model[i][j] == JKR) use_history = 1;



  size_history = 3*tangential_history + 3*roll_history + twist_history;
@@ -920,6 +942,19 @@ void PairGranular::init_style() 
      else twist_history_index = 0;

    }

  }

  for (int i = 1; i <= atom->ntypes; i++)

    for (int j = i; j <= atom->ntypes; j++)

      if (tangential_model[i][j] == TANGENTIAL_MINDLIN_RESCALE){

        size_history += 1;

        roll_history_index += 1;

        twist_history_index += 1;

        nondefault_history_transfer = 1;

        history_transfer_factors = new int[size_history];

        for (int ii = 0; ii < size_history; ++ii)

          history_transfer_factors[ii] = -1;

        history_transfer_factors[3] = 1;

        break;

      }



  int irequest = neighbor->request(this,instance_me);

  neighbor->requests[irequest]->size = 1;
@@ -1612,3 +1647,12 @@ double PairGranular::pulloff_distance(double radi, double radj, int itype, int j 
  return a*a/Reff - 2*sqrt(M_PI*coh*a/E);

}



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

     Transfer history during fix/neigh/history exchange

      Only needed if any history entries i-j are not just negative of j-i entries

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

void PairGranular::transfer_history(double* source, double* target){

  for (int i = 0; i < size_history; i++)

    target[i] = history_transfer_factors[i]*source[i];

}

src/GRANULAR/pair_granular.h

+2 −0
Original line number Diff line number Diff line
@@ -61,9 +61,11 @@ public: 
  int nmax;                // allocated size of mass_rigid



  void allocate();

  void transfer_history(double*, double*);



private:

  int size_history;

  int *history_transfer_factors;



  //Model choices

  int **normal_model, **damping_model;

CRA Git | Maintained and supported by SUSTech CRA and CCSE