Commit 921b6d81 authored by Jibril B. Coulibaly's avatar Jibril B. Coulibaly
Browse files

relative threshold for contact frame update based on tangential critical force

parent ee6ef98b

src/GRANULAR/pair_granular.cpp

+48 −38
Original line number Diff line number Diff line
@@ -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,6 +442,13 @@ 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 (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
@@ -456,7 +464,7 @@ void PairGranular::compute(int eflag, int vflag) 
              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,9 +561,14 @@ 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;

          k_roll = roll_coeffs[itype][jtype][0];

          damp_roll = roll_coeffs[itype][jtype][1];

          Frcrit = roll_coeffs[itype][jtype][2] * Fncrit;



          if (historyupdate) {

            // rotate into tangential plane

            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]);
@@ -572,20 +585,17 @@ 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;

          }



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

CRA Git | Maintained and supported by SUSTech CRA and CCSE