diff --git a/doc/src/pair_granular.txt b/doc/src/pair_granular.txt index ff2ee94be0..911e3cc1dc 100644 --- a/doc/src/pair_granular.txt +++ b/doc/src/pair_granular.txt @@ -533,14 +533,16 @@ the pairwise interaction force. However, the single() function also calculates 10 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 relative velocity in the -normal direction (along the line joining the 2 sphere centers). The -last 3 (8-10) the components of the relative velocity in the -tangential direction. +The next 3 (5-7) are the components of the rolling torque acting on +particle I. The next entry (8) is the magnitude of the rolling torque. +The next entry (9) is the magnitude of the twisting torque acting +about the vector connecting the two particle centers. +The last 3 (10-12) are the components of the vector connecting +the centers of the two particles (x_I - x_J). These extra quantities can be accessed by the "compute pair/local"_compute_pair_local.html command, as {p1}, {p2}, ..., -{p10}. +{p12}. :line @@ -568,6 +570,7 @@ compute depend on atom velocities. See the [Related commands:] "pair_coeff"_pair_coeff.html +"pair gran/*"_pair_gran.html [Default:] diff --git a/src/GRANULAR/fix_wall_gran.cpp b/src/GRANULAR/fix_wall_gran.cpp index d3129b7cdb..6e8cba7b4f 100644 --- a/src/GRANULAR/fix_wall_gran.cpp +++ b/src/GRANULAR/fix_wall_gran.cpp @@ -39,11 +39,24 @@ using namespace MathConst; // XYZ PLANE need to be 0,1,2 enum{XPLANE=0,YPLANE=1,ZPLANE=2,ZCYLINDER,REGION}; -enum{HOOKE,HOOKE_HISTORY,HERTZ_HISTORY,JKR_ROLLING,DMT_ROLLING}; +enum{HOOKE,HOOKE_HISTORY,HERTZ_HISTORY,GRANULAR}; enum{NONE,CONSTANT,EQUAL}; -enum {TSUJI, BRILLIANTOV}; -enum {INDEP, BRILLROLL}; +#define PI27SQ 266.47931882941264802866 // 27*PI**2 +#define THREEROOT3 5.19615242270663202362 // 3*sqrt(3) +#define SIXROOT6 14.69693845669906728801 // 6*sqrt(6) +#define INVROOT6 0.40824829046386307274 // 1/sqrt(6) +#define FOURTHIRDS 1.333333333333333 // 4/3 +#define THREEQUARTERS 0.75 // 3/4 +#define TWOPI 6.28318530717959 // 2*PI + +#define EPSILON 1e-10 + +enum {NORMAL_HOOKE, NORMAL_HERTZ, HERTZ_MATERIAL, DMT, JKR}; +enum {VELOCITY, VISCOELASTIC, TSUJI}; +enum {TANGENTIAL_NOHISTORY, TANGENTIAL_HISTORY, TANGENTIAL_MINDLIN}; +enum {TWIST_NONE, TWIST_SDS, TWIST_MARSHALL}; +enum {ROLL_NONE, ROLL_SDS}; #define BIG 1.0e20 #define EPSILON 1e-10 @@ -51,7 +64,7 @@ enum {INDEP, BRILLROLL}; /* ---------------------------------------------------------------------- */ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), idregion(NULL), shearone(NULL), fix_rigid(NULL), mass_rigid(NULL) + Fix(lmp, narg, arg), idregion(NULL), history_one(NULL), fix_rigid(NULL), mass_rigid(NULL) { if (narg < 4) error->all(FLERR,"Illegal fix wall/gran command"); @@ -66,19 +79,18 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : if (strcmp(arg[3],"hooke") == 0) pairstyle = HOOKE; else if (strcmp(arg[3],"hooke/history") == 0) pairstyle = HOOKE_HISTORY; else if (strcmp(arg[3],"hertz/history") == 0) pairstyle = HERTZ_HISTORY; - else if (strcmp(arg[3],"dmt/rolling") == 0) pairstyle = DMT_ROLLING; - //else if (strcmp(arg[3],"jkr/rolling") == 0) pairstyle = JKR_ROLLING; + else if (strcmp(arg[3],"granular") == 0) pairstyle = GRANULAR; else error->all(FLERR,"Invalid fix wall/gran interaction style"); - history = restart_peratom = 1; - if (pairstyle == HOOKE) history = restart_peratom = 0; + use_history = restart_peratom = 1; + if (pairstyle == HOOKE) use_history = restart_peratom = 0; // wall/particle coefficients int iarg; - if (pairstyle != JKR_ROLLING && pairstyle != DMT_ROLLING) { - sheardim = 3; + if (pairstyle != GRANULAR) { + size_history = 3; if (narg < 11) error->all(FLERR,"Illegal fix wall/gran command"); kn = force->numeric(FLERR,arg[4]); @@ -103,46 +115,157 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : kn /= force->nktv2p; kt /= force->nktv2p; } - iarg = 10; } else { - if (narg < 12) error->all(FLERR,"Illegal fix wall/gran command"); - - sheardim = 7; - Emod = force->numeric(FLERR,arg[4]); - Gmod = force->numeric(FLERR,arg[5]); - xmu = force->numeric(FLERR,arg[6]); - gamman = force->numeric(FLERR,arg[7]); - Ecoh = force->numeric(FLERR,arg[8]); - kR = force->numeric(FLERR,arg[9]); - muR = force->numeric(FLERR,arg[10]); - etaR = force->numeric(FLERR,arg[11]); - - //Defaults - normaldamp = TSUJI; - rollingdamp = INDEP; - - iarg = 12; - for (int iiarg=iarg; iiarg < narg; ++iiarg){ - if (strcmp(arg[iiarg], "normaldamp") == 0){ - if(iiarg+2 > narg) error->all(FLERR, "Invalid fix/wall/gran region command"); - if (strcmp(arg[iiarg+1],"tsuji") == 0){ - normaldamp = TSUJI; - alpha = gamman; - } - else if (strcmp(arg[iiarg+1],"brilliantov") == 0) normaldamp = BRILLIANTOV; - else error->all(FLERR, "Invalid normal damping model for fix wall/gran dmt/rolling"); - iarg += 2; + iarg = 4; + damping_model = VISCOELASTIC; + roll_model = twist_model = NONE; + while (iarg < narg){ + if (strcmp(arg[iarg], "hooke") == 0){ + if (iarg + 2 >= narg) error->all(FLERR,"Illegal fix wall/gran command, not enough parameters provided for Hooke option"); + normal_model = NORMAL_HOOKE; + normal_coeffs[0] = force->numeric(FLERR,arg[iarg+1]); //kn + normal_coeffs[1] = force->numeric(FLERR,arg[iarg+2]); //damping + iarg += 3; } - if (strcmp(arg[iiarg], "rollingdamp") == 0){ - if(iiarg+2 > narg) error->all(FLERR, "Invalid fix/wall/gran region command"); - if (strcmp(arg[iarg+1],"independent") == 0) rollingdamp = INDEP; - else if (strcmp(arg[iarg+1],"brilliantov") == 0) rollingdamp = BRILLROLL; - else error->all(FLERR, "Invalid rolling damping model for fix wall/gran dmt/rolling"); - iarg += 2; + else if (strcmp(arg[iarg], "hertz") == 0){ + int num_coeffs = 2; + if (iarg + num_coeffs >= narg) error->all(FLERR,"Illegal fix wall/gran command, not enough parameters provided for Hertz option"); + normal_model = NORMAL_HERTZ; + normal_coeffs[0] = force->numeric(FLERR,arg[iarg+1]); //kn + normal_coeffs[1] = force->numeric(FLERR,arg[iarg+2]); //damping + iarg += num_coeffs+1; + } + else if (strcmp(arg[iarg], "hertz/material") == 0){ + int num_coeffs = 3; + if (iarg + num_coeffs >= narg) error->all(FLERR,"Illegal fix wall/gran command, not enough parameters provided for Hertz option"); + normal_model = HERTZ_MATERIAL; + Emod = force->numeric(FLERR,arg[iarg+1]); //E + normal_coeffs[1] = force->numeric(FLERR,arg[iarg+2]); //damping + poiss = force->numeric(FLERR,arg[iarg+3]); //Poisson's ratio + normal_coeffs[0] = Emod/(2*(1-poiss))*FOURTHIRDS; + normal_coeffs[2] = poiss; + iarg += num_coeffs+1; + } + else if (strcmp(arg[iarg], "dmt") == 0){ + if (iarg + 4 >= narg) error->all(FLERR,"Illegal fix wall/gran command, not enough parameters provided for Hertz option"); + normal_model = DMT; + Emod = force->numeric(FLERR,arg[iarg+1]); //E + normal_coeffs[1] = force->numeric(FLERR,arg[iarg+2]); //damping + poiss = force->numeric(FLERR,arg[iarg+3]); //Poisson's ratio + normal_coeffs[0] = Emod/(2*(1-poiss))*FOURTHIRDS; + normal_coeffs[2] = poiss; + normal_coeffs[3] = force->numeric(FLERR,arg[iarg+4]); //cohesion + iarg += 5; + } + else if (strcmp(arg[iarg], "jkr") == 0){ + if (iarg + 4 >= narg) error->all(FLERR,"Illegal wall/gran command, not enough parameters provided for JKR option"); + beyond_contact = 1; + normal_model = JKR; + Emod = force->numeric(FLERR,arg[iarg+1]); //E + normal_coeffs[1] = force->numeric(FLERR,arg[iarg+2]); //damping + poiss = force->numeric(FLERR,arg[iarg+3]); //Poisson's ratio + normal_coeffs[0] = Emod/(2*(1-poiss))*FOURTHIRDS; + normal_coeffs[2] = poiss; + normal_coeffs[3] = force->numeric(FLERR,arg[iarg+4]); //cohesion + iarg += 5; + } + else if (strcmp(arg[iarg], "damping") == 0){ + if (iarg+1 >= narg) error->all(FLERR, "Illegal wall/gran command, not enough parameters provided for damping model"); + if (strcmp(arg[iarg+1], "velocity") == 0){ + damping_model = VELOCITY; + iarg += 1; + } + else if (strcmp(arg[iarg+1], "viscoelastic") == 0){ + damping_model = VISCOELASTIC; + iarg += 1; + } + else if (strcmp(arg[iarg+1], "tsuji") == 0){ + damping_model = TSUJI; + iarg += 1; + } + else error->all(FLERR, "Illegal wall/gran command, unrecognized damping model"); + iarg += 1; + } + else if (strcmp(arg[iarg], "tangential") == 0){ + if (iarg + 1 >= narg) error->all(FLERR,"Illegal pair_coeff command, must specify tangential model after 'tangential' keyword"); + if (strcmp(arg[iarg+1], "linear_nohistory") == 0){ + if (iarg + 3 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for tangential model"); + tangential_model = TANGENTIAL_NOHISTORY; + tangential_coeffs[0] = 0; + tangential_coeffs[1] = force->numeric(FLERR,arg[iarg+2]); //gammat + tangential_coeffs[2] = force->numeric(FLERR,arg[iarg+3]); //friction coeff. + iarg += 4; + } + else if (strcmp(arg[iarg+1], "linear_history") == 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; + tangential_coeffs[0] = force->numeric(FLERR,arg[iarg+2]); //kt + tangential_coeffs[1] = force->numeric(FLERR,arg[iarg+3]); //gammat + tangential_coeffs[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff. + iarg += 5; + } + else{ + error->all(FLERR, "Illegal pair_coeff command, tangential model not recognized"); + } + } + else if (strcmp(arg[iarg], "rolling") == 0){ + if (iarg + 1 >= narg) error->all(FLERR, "Illegal wall/gran command, not enough parameters"); + if (strcmp(arg[iarg+1], "none") == 0){ + roll_model = ROLL_NONE; + iarg += 2; + } + else if (strcmp(arg[iarg+1], "sds") == 0){ + if (iarg + 4 >= narg) error->all(FLERR,"Illegal wall/gran command, not enough parameters provided for rolling model"); + roll_model = ROLL_SDS; + roll_history = 1; + roll_coeffs[0] = force->numeric(FLERR,arg[iarg+2]); //kR + roll_coeffs[1] = force->numeric(FLERR,arg[iarg+3]); //gammaR + roll_coeffs[2] = force->numeric(FLERR,arg[iarg+4]); //rolling friction coeff. + iarg += 5; + } + else{ + error->all(FLERR, "Illegal wall/gran command, rolling friction model not recognized"); + } + } + else if (strcmp(arg[iarg], "twisting") == 0){ + if (iarg + 1 >= narg) error->all(FLERR, "Illegal wall/gran command, not enough parameters"); + if (strcmp(arg[iarg+1], "none") == 0){ + twist_model = TWIST_NONE; + iarg += 2; + } + else if (strcmp(arg[iarg+1], "marshall") == 0){ + twist_model = TWIST_MARSHALL; + twist_history = 1; + iarg += 2; + } + else if (strcmp(arg[iarg+1], "sds") == 0){ + if (iarg + 4 >= narg) error->all(FLERR,"Illegal wall/gran command, not enough parameters provided for twist model"); + twist_model = TWIST_SDS; + twist_history = 1; + twist_coeffs[0] = force->numeric(FLERR,arg[iarg+2]); //kt + twist_coeffs[1] = force->numeric(FLERR,arg[iarg+3]); //gammat + twist_coeffs[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff. + iarg += 5; + } + else{ + error->all(FLERR, "Illegal wall/gran command, twisting friction model not recognized"); + } + } + else if (strcmp(arg[iarg], "xplane") == 0 || + strcmp(arg[iarg], "yplane") == 0 || + strcmp(arg[iarg], "zplane") == 0 || + strcmp(arg[iarg], "zcylinder") == 0 || + strcmp(arg[iarg], "region") == 0){ + break; + } + else{ + error->all(FLERR, "Illegal fix wall/gran command"); } } + size_history = (normal_model == JKR) + 3*tangential_history + 3*roll_history + twist_history; } // wallstyle args @@ -216,11 +339,10 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : iarg += 3; } else if (strcmp(arg[iarg],"store_contacts") == 0){ peratom_flag = 1; - size_peratom_cols = 8; //Could make this a user input option? + size_peratom_cols = 8; peratom_freq = 1; iarg += 1; } else error->all(FLERR,"Illegal fix wall/gran command"); - } if (wallstyle == XPLANE && domain->xperiodic) @@ -252,7 +374,7 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : // perform initial allocation of atom-based arrays // register with Atom class - shearone = NULL; + history_one = NULL; grow_arrays(atom->nmax); atom->add_callback(0); atom->add_callback(1); @@ -260,14 +382,14 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : nmax = 0; mass_rigid = NULL; - // initialize shear history as if particle is not touching region - // shearone will be NULL for wallstyle = REGION + // initialize history as if particle is not touching region + // history_one will be NULL for wallstyle = REGION - if (history && shearone) { + if (use_history && history_one) { int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) - for (int j = 0; j < sheardim; j++) - shearone[i][j] = 0.0; + for (int j = 0; j < size_history; j++) + history_one[i][j] = 0.0; } if (peratom_flag){ @@ -292,7 +414,7 @@ FixWallGran::~FixWallGran() // delete local storage delete [] idregion; - memory->destroy(shearone); + memory->destroy(history_one); memory->destroy(mass_rigid); } @@ -323,6 +445,34 @@ void FixWallGran::init() for (i = 0; i < modify->nfix; i++) if (modify->fix[i]->rigid_flag) break; if (i < modify->nfix) fix_rigid = modify->fix[i]; + + tangential_history_index = 0; + if (roll_history){ + if (tangential_history) roll_history_index = 3; + else roll_history_index = 0; + } + if (twist_history){ + if (tangential_history){ + if (roll_history) twist_history_index = 6; + else twist_history_index = 3; + } + else{ + if (roll_history) twist_history_index = 3; + else twist_history_index = 0; + } + } + if (normal_model == JKR){ + tangential_history_index += 1; + roll_history_index += 1; + twist_history_index += 1; + } + + if (damping_model == TSUJI){ + double cor = normal_coeffs[1]; + normal_coeffs[1] = 1.2728-4.2783*cor+11.087*pow(cor,2)-22.348*pow(cor,3)+ + 27.467*pow(cor,4)-18.022*pow(cor,5)+ + 4.8218*pow(cor,6); + } } /* ---------------------------------------------------------------------- */ @@ -346,10 +496,10 @@ void FixWallGran::post_force(int /*vflag*/) double dx,dy,dz,del1,del2,delxy,delr,rsq,rwall,meff; double vwall[3]; - // do not update shear history during setup + // do not update history during setup - shearupdate = 1; - if (update->setupflag) shearupdate = 0; + history_update = 1; + if (update->setupflag) history_update = 0; // if just reneighbored: // update rigid body masses for owned atoms if using FixRigid @@ -395,7 +545,7 @@ void FixWallGran::post_force(int /*vflag*/) // if wall was set to NULL, it's skipped since lo/hi are infinity // compute force and torque on atom if close enough to wall // via wall potential matched to pair potential - // set shear if pair potential stores history + // set history if pair potential stores history double **x = atom->x; double **v = atom->v; @@ -450,12 +600,27 @@ void FixWallGran::post_force(int /*vflag*/) rsq = dx*dx + dy*dy + dz*dz; - if (rsq > radius[i]*radius[i]) { - if (history) - for (j = 0; j < sheardim; j++) - shearone[i][j] = 0.0; + double rad; + if (pairstyle == GRANULAR && normal_model == JKR){ + rad = radius[i] + pulloff_distance(radius[i]); + } + else + rad = radius[i]; - } else { + if (rsq > rad*rad) { + if (use_history) + for (j = 0; j < size_history; j++) + history_one[i][j] = 0.0; + } + else { + if (pairstyle == GRANULAR && normal_model == JKR && use_history){ + if ((history_one[i][0] == 0) && (rsq > radius[i]*radius[i])){ + // Particles have not contacted yet, and are outside of contact distance + for (j = 0; j < size_history; j++) + history_one[i][j] = 0.0; + continue; + } + } // meff = effective mass of sphere // if I is part of rigid body, use body mass @@ -484,20 +649,16 @@ void FixWallGran::post_force(int /*vflag*/) omega[i],torque[i],radius[i],meff, contact); else if (pairstyle == HOOKE_HISTORY) hooke_history(rsq,dx,dy,dz,vwall,v[i],f[i], - omega[i],torque[i],radius[i],meff,shearone[i], + omega[i],torque[i],radius[i],meff,history_one[i], contact); else if (pairstyle == HERTZ_HISTORY) hertz_history(rsq,dx,dy,dz,vwall,rwall,v[i],f[i], - omega[i],torque[i],radius[i],meff,shearone[i], + omega[i],torque[i],radius[i],meff,history_one[i], contact); - else if (pairstyle == DMT_ROLLING) - dmt_rolling(rsq,dx,dy,dz,vwall,rwall,v[i],f[i], - omega[i],torque[i],radius[i],meff,shearone[i], + else if (pairstyle == GRANULAR) + granular(rsq,dx,dy,dz,vwall,rwall,v[i],f[i], + omega[i],torque[i],radius[i],meff,history_one[i], contact); - /*else if (pairstyle == JKR_ROLLING) - jkr_rolling(rsq,dx,dy,dz,vwall,rwall,v[i],f[i], - omega[i],torque[i],radius[i],meff,shearone[i], - contact);*/ } } } @@ -605,7 +766,7 @@ void FixWallGran::hooke(double rsq, double dx, double dy, double dz, void FixWallGran::hooke_history(double rsq, double dx, double dy, double dz, double *vwall, double *v, double *f, double *omega, double *torque, - double radius, double meff, double *shear, + double radius, double meff, double *history, double *contact) { double r,vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; @@ -657,28 +818,28 @@ void FixWallGran::hooke_history(double rsq, double dx, double dy, double dz, // shear history effects - if (shearupdate) { - shear[0] += vtr1*dt; - shear[1] += vtr2*dt; - shear[2] += vtr3*dt; + if (history_update) { + history[0] += vtr1*dt; + history[1] += vtr2*dt; + history[2] += vtr3*dt; } - shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + shear[2]*shear[2]); + shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + history[2]*history[2]); // rotate shear displacements - rsht = shear[0]*dx + shear[1]*dy + shear[2]*dz; + rsht = history[0]*dx + history[1]*dy + history[2]*dz; rsht = rsht*rsqinv; - if (shearupdate) { - shear[0] -= rsht*dx; - shear[1] -= rsht*dy; - shear[2] -= rsht*dz; + if (history_update) { + history[0] -= rsht*dx; + history[1] -= rsht*dy; + history[2] -= rsht*dz; } // tangential forces = shear + tangential velocity damping - fs1 = - (kt*shear[0] + meff*gammat*vtr1); - fs2 = - (kt*shear[1] + meff*gammat*vtr2); - fs3 = - (kt*shear[2] + meff*gammat*vtr3); + fs1 = - (kt*history[0] + meff*gammat*vtr1); + fs2 = - (kt*history[1] + meff*gammat*vtr2); + fs3 = - (kt*history[2] + meff*gammat*vtr3); // rescale frictional displacements and forces if needed @@ -687,11 +848,11 @@ void FixWallGran::hooke_history(double rsq, double dx, double dy, double dz, if (fs > fn) { if (shrmag != 0.0) { - shear[0] = (fn/fs) * (shear[0] + meff*gammat*vtr1/kt) - + history[0] = (fn/fs) * (history[0] + meff*gammat*vtr1/kt) - meff*gammat*vtr1/kt; - shear[1] = (fn/fs) * (shear[1] + meff*gammat*vtr2/kt) - + history[1] = (fn/fs) * (history[1] + meff*gammat*vtr2/kt) - meff*gammat*vtr2/kt; - shear[2] = (fn/fs) * (shear[2] + meff*gammat*vtr3/kt) - + history[2] = (fn/fs) * (history[2] + meff*gammat*vtr3/kt) - meff*gammat*vtr3/kt; fs1 *= fn/fs ; fs2 *= fn/fs; @@ -728,7 +889,7 @@ void FixWallGran::hooke_history(double rsq, double dx, double dy, double dz, void FixWallGran::hertz_history(double rsq, double dx, double dy, double dz, double *vwall, double rwall, double *v, double *f, double *omega, double *torque, - double radius, double meff, double *shear, + double radius, double meff, double *history, double *contact) { double r,vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; @@ -787,28 +948,28 @@ void FixWallGran::hertz_history(double rsq, double dx, double dy, double dz, // shear history effects - if (shearupdate) { - shear[0] += vtr1*dt; - shear[1] += vtr2*dt; - shear[2] += vtr3*dt; + if (history_update) { + history[0] += vtr1*dt; + history[1] += vtr2*dt; + history[2] += vtr3*dt; } - shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + shear[2]*shear[2]); + shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + history[2]*history[2]); - // rotate shear displacements + // rotate history displacements - rsht = shear[0]*dx + shear[1]*dy + shear[2]*dz; + rsht = history[0]*dx + history[1]*dy + history[2]*dz; rsht = rsht*rsqinv; - if (shearupdate) { - shear[0] -= rsht*dx; - shear[1] -= rsht*dy; - shear[2] -= rsht*dz; + if (history_update) { + history[0] -= rsht*dx; + history[1] -= rsht*dy; + history[2] -= rsht*dz; } // tangential forces = shear + tangential velocity damping - fs1 = -polyhertz * (kt*shear[0] + meff*gammat*vtr1); - fs2 = -polyhertz * (kt*shear[1] + meff*gammat*vtr2); - fs3 = -polyhertz * (kt*shear[2] + meff*gammat*vtr3); + fs1 = -polyhertz * (kt*history[0] + meff*gammat*vtr1); + fs2 = -polyhertz * (kt*history[1] + meff*gammat*vtr2); + fs3 = -polyhertz * (kt*history[2] + meff*gammat*vtr3); // rescale frictional displacements and forces if needed @@ -817,11 +978,11 @@ void FixWallGran::hertz_history(double rsq, double dx, double dy, double dz, if (fs > fn) { if (shrmag != 0.0) { - shear[0] = (fn/fs) * (shear[0] + meff*gammat*vtr1/kt) - + history[0] = (fn/fs) * (history[0] + meff*gammat*vtr1/kt) - meff*gammat*vtr1/kt; - shear[1] = (fn/fs) * (shear[1] + meff*gammat*vtr2/kt) - + history[1] = (fn/fs) * (history[1] + meff*gammat*vtr2/kt) - meff*gammat*vtr2/kt; - shear[2] = (fn/fs) * (shear[2] + meff*gammat*vtr3/kt) - + history[2] = (fn/fs) * (history[2] + meff*gammat*vtr3/kt) - meff*gammat*vtr3/kt; fs1 *= fn/fs ; fs2 *= fn/fs; @@ -854,268 +1015,366 @@ void FixWallGran::hertz_history(double rsq, double dx, double dy, double dz, } -void FixWallGran::dmt_rolling(double rsq, double dx, double dy, double dz, +void FixWallGran::granular(double rsq, double dx, double dy, double dz, double *vwall, double rwall, double *v, double *f, double *omega, double *torque, - double radius, double meff, double *shear, + double radius, double meff, double *history, double *contact) { - int i,j,ii,jj,inum,jnum; - int itype,jtype; - double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz,nx,ny,nz; - double radi,radj,radsum,r,rinv,rsqinv,R,a; - double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; - double wr1,wr2,wr3; - double vtr1,vtr2,vtr3,vrel; - double kn, kt, k_Q, k_R, eta_N, eta_T, eta_Q, eta_R; - double Fhz, Fdamp, Fdmt, Fne, Fntot, Fscrit, Frcrit; - double overlap; - double mi,mj,damp,ccel,tor1,tor2,tor3; - double relrot1,relrot2,relrot3,vrl1,vrl2,vrl3,vrlmag,vrlmaginv; - double rollmag, rolldotn, scalefac; - double fr, fr1, fr2, fr3; - double signtwist, magtwist, magtortwist, Mtcrit; - double fs,fs1,fs2,fs3,roll1,roll2,roll3,torroll1,torroll2,torroll3; - double tortwist1, tortwist2, tortwist3; - double shrmag,rsht; + 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,r,rinv,rsqinv; + double Reff, delta, dR, dR2; - r = sqrt(rsq); - rinv = 1.0/r; - rsqinv = 1.0/rsq; + double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; + double wr1,wr2,wr3; + double vtr1,vtr2,vtr3,vrel; - radsum = radius + rwall; - if (rwall == 0) R = radius; - else R = radius*rwall/(radius+rwall); + double knfac, damp_normal, damp_normal_prefactor; + double k_tangential, damp_tangential; + double Fne, Ft, Fdamp, Fntot, Fncrit, Fscrit, Frcrit; + double fs, fs1, fs2, fs3; - nx = dx*rinv; - ny = dy*rinv; - nz = dz*rinv; + double mi,mj,damp,ccel,tor1,tor2,tor3; + double relrot1,relrot2,relrot3,vrl1,vrl2,vrl3,vrlmag,vrlmaginv; - // relative translational velocity + //For JKR + double R2, coh, F_pulloff, delta_pulloff, dist_pulloff, a, a2, E; + double t0, t1, t2, t3, t4, t5, t6; + double sqrt1, sqrt2, sqrt3, sqrt4; - vr1 = v[0] - vwall[0]; - vr2 = v[1] - vwall[1]; - vr3 = v[2] - vwall[2]; + //Rolling + double k_roll, damp_roll; + double roll1, roll2, roll3, torroll1, torroll2, torroll3; + double rollmag, rolldotn, scalefac; + double fr, fr1, fr2, fr3; - // normal component + //Twisting + double k_twist, damp_twist, mu_twist; + double signtwist, magtwist, magtortwist, Mtcrit; + double tortwist1, tortwist2, tortwist3; - vnnr = vr1*nx + vr2*ny + vr3*nz; //v_R . n - vn1 = nx*vnnr; - vn2 = ny*vnnr; - vn3 = nz*vnnr; + double shrmag,rsht; + int *ilist,*jlist,*numneigh,**firstneigh; + int *touch,**firsttouch; + double *allhistory,**firsthistory; + r = sqrt(rsq); + radsum = rwall + radius; - //**************************************** - //Normal force = Hertzian contact + DMT + damping - //**************************************** - overlap = radsum - r; - a = sqrt(R*overlap); - kn = 4.0/3.0*Emod*a; - Fhz = kn*overlap; + E = normal_coeffs[0]; - //Damping (based on Tsuji et al) - if (normaldamp == BRILLIANTOV) eta_N = a*meff*gamman; - else if (normaldamp == TSUJI) eta_N=alpha*sqrt(meff*kn); + radsum = radius + rwall; + if (rwall == 0) Reff = radius; + else Reff = radius*rwall/(radius+rwall); - Fdamp = -eta_N*vnnr; //F_nd eq 23 and Zhao eq 19 + rinv = 1.0/r; - //DMT - Fdmt = -4*MY_PI*Ecoh*R; + nx = dx*rinv; + ny = dy*rinv; + nz = dz*rinv; - Fne = Fhz + Fdmt; - Fntot = Fne + Fdamp; + // relative translational velocity - //**************************************** - //Tangential force, including shear history effects - //**************************************** + vr1 = v[0] - vwall[0]; + vr2 = v[1] - vwall[1]; + vr3 = v[2] - vwall[2]; - // tangential component - vt1 = vr1 - vn1; - vt2 = vr2 - vn2; - vt3 = vr3 - vn3; + // normal component - // relative rotational velocity + vnnr = vr1*nx + vr2*ny + vr3*nz; //v_R . n + vn1 = nx*vnnr; + vn2 = ny*vnnr; + vn3 = nz*vnnr; - wr1 = radius*omega[0] * rinv; - wr2 = radius*omega[1] * rinv; - wr3 = radius*omega[2] * rinv; + delta = radsum - r; + dR = delta*Reff; + if (normal_model == JKR){ + history[0] = 1.0; + E *= THREEQUARTERS; + R2=Reff*Reff; + coh = normal_coeffs[3]; + dR2 = dR*dR; + t0 = coh*coh*R2*R2*E; + t1 = PI27SQ*t0; + t2 = 8*dR*dR2*E*E*E; + t3 = 4*dR2*E; + sqrt1 = MAX(0, t0*(t1+2*t2)); //In case of sqrt(0) < 0 due to precision issues + t4 = cbrt(t1+t2+THREEROOT3*M_PI*sqrt(sqrt1)); + t5 = t3/t4 + t4/E; + sqrt2 = MAX(0, 2*dR + t5); + t6 = sqrt(sqrt2); + sqrt3 = MAX(0, 4*dR - t5 + SIXROOT6*coh*M_PI*R2/(E*t6)); + a = INVROOT6*(t6 + sqrt(sqrt3)); + a2 = a*a; + knfac = normal_coeffs[0]*a; + Fne = knfac*a2/Reff - TWOPI*a2*sqrt(4*coh*E/(M_PI*a)); + } + else{ + knfac = E; //Hooke + a = sqrt(dR); + if (normal_model != HOOKE){ + Fne *= a; + knfac *= a; + } + Fne = knfac*delta; + if (normal_model == DMT) + Fne -= 4*MY_PI*normal_coeffs[3]*Reff; + } - // relative tangential velocities - vtr1 = vt1 - (nz*wr2-ny*wr3); - vtr2 = vt2 - (nx*wr3-nz*wr1); - vtr3 = vt3 - (ny*wr1-nx*wr2); - vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3; - vrel = sqrt(vrel); + if (damping_model == VELOCITY){ + damp_normal = 1; + } + else if (damping_model == VISCOELASTIC){ + damp_normal = a*meff; + } + else if (damping_model == TSUJI){ + damp_normal = sqrt(meff*knfac); + } - // shear history effects - shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + - shear[2]*shear[2]); + damp_normal_prefactor = normal_coeffs[1]*damp_normal; + Fdamp = -damp_normal_prefactor*vnnr; - // Rotate and update shear displacements. - // See e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235 - if (shearupdate) { - rsht = shear[0]*nx + shear[1]*ny + shear[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! - shear[0] -= rsht*nx; - shear[1] -= rsht*ny; - shear[2] -= rsht*nz; - //Also rescale to preserve magnitude - shear[0] *= scalefac; - shear[1] *= scalefac; - shear[2] *= scalefac; - } - //Update shear history - shear[0] += vtr1*dt; - shear[1] += vtr2*dt; - shear[2] += vtr3*dt; - } + Fntot = Fne + Fdamp; - // tangential forces = shear + tangential velocity damping - // following Zhao and Marshall Phys Fluids v20, p043302 (2008) - kt=8.0*Gmod*a; + //**************************************** + //Tangential force, including history effects + //**************************************** - eta_T = eta_N; //Based on discussion in Marshall; eta_T can also be an independent parameter - fs1 = -kt*shear[0] - eta_T*vtr1; //eq 26 - fs2 = -kt*shear[1] - eta_T*vtr2; - fs3 = -kt*shear[2] - eta_T*vtr3; + // tangential component + vt1 = vr1 - vn1; + vt2 = vr2 - vn2; + vt3 = vr3 - vn3; - // rescale frictional displacements and forces if needed - Fscrit = xmu * fabs(Fne); - // For JKR, use eq 43 of Marshall. For DMT, use Fne instead + // relative rotational velocity + wr1 = radius*omega[0] * rinv; + wr2 = radius*omega[1] * rinv; + wr3 = radius*omega[2] * rinv; - shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + shear[2]*shear[2]); + // relative tangential velocities + vtr1 = vt1 - (nz*wr2-ny*wr3); + vtr2 = vt2 - (nx*wr3-nz*wr1); + vtr3 = vt3 - (ny*wr1-nx*wr2); + vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3; + vrel = sqrt(vrel); - fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); - if (fs > Fscrit) { - if (shrmag != 0.0) { - //shear[0] = (Fcrit/fs) * (shear[0] + eta_T*vtr1/kt) - eta_T*vtr1/kt; - //shear[1] = (Fcrit/fs) * (shear[1] + eta_T*vtr1/kt) - eta_T*vtr1/kt; - //shear[2] = (Fcrit/fs) * (shear[2] + eta_T*vtr1/kt) - eta_T*vtr1/kt; - shear[0] = -1.0/kt*(Fscrit*fs1/fs + eta_T*vtr1); //Same as above, but simpler (check!) - shear[1] = -1.0/kt*(Fscrit*fs2/fs + eta_T*vtr2); - shear[2] = -1.0/kt*(Fscrit*fs3/fs + eta_T*vtr3); - fs1 *= Fscrit/fs; - fs2 *= Fscrit/fs; - fs3 *= Fscrit/fs; - } else fs1 = fs2 = fs3 = 0.0; - } + if (normal_model == JKR){ + F_pulloff = 3*M_PI*coh*Reff; + Fncrit = fabs(Fne + 2*F_pulloff); + } + else if (normal_model == DMT){ + F_pulloff = 4*M_PI*coh*Reff; + Fncrit = fabs(Fne + 2*F_pulloff); + } + else{ + Fncrit = fabs(Fntot); + } - //**************************************** - // Rolling force, including shear history effects - //**************************************** + //------------------------------ + //Tangential forces + //------------------------------ + k_tangential = tangential_coeffs[0]; + damp_tangential = tangential_coeffs[1]*damp_normal_prefactor; - relrot1 = omega[0]; //- omega[j][0]; TODO: figure out how to - relrot2 = omega[1]; //- omega[j][1]; incorporate wall angular - relrot3 = omega[2]; //- omega[j][2]; velocity + int thist0 = tangential_history_index; + int thist1 = thist0 + 1; + int thist2 = thist1 + 1; - // 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 = R*(relrot2*nz - relrot3*ny); //- 0.5*((radj-radi)/radsum)*vtr1; - vrl2 = R*(relrot3*nx - relrot1*nz); //- 0.5*((radj-radi)/radsum)*vtr2; - vrl3 = R*(relrot1*ny - relrot2*nx); //- 0.5*((radj-radi)/radsum)*vtr3; - vrlmag = sqrt(vrl1*vrl1+vrl2*vrl2+vrl3*vrl3); - if (vrlmag != 0.0) vrlmaginv = 1.0/vrlmag; - else vrlmaginv = 0.0; + if (tangential_history){ + shrmag = sqrt(history[thist0]*history[thist0] + history[thist1]*history[thist1] + + history[thist2]*history[thist2]); - // Rolling displacement - rollmag = sqrt(shear[3]*shear[3] + shear[4]*shear[4] + shear[5]*shear[5]); - rolldotn = shear[3]*nx + shear[4]*ny + shear[5]*nz; + // 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 (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! + history[thist0] -= rsht*nx; + history[thist1] -= rsht*ny; + history[thist2] -= rsht*nz; + //Also rescale to preserve magnitude + history[thist0] *= scalefac; + history[thist1] *= scalefac; + history[thist2] *= scalefac; + } + //Update history + history[thist0] += vtr1*dt; + history[thist1] += vtr2*dt; + history[thist2] += vtr3*dt; + } - if (shearupdate) { - if (fabs(rolldotn) < EPSILON) rolldotn = 0; - if (rolldotn > 0){ //Rotate into tangential plane - scalefac = rollmag/(rollmag - rolldotn); - shear[3] -= rolldotn*nx; - shear[4] -= rolldotn*ny; - shear[5] -= rolldotn*nz; - //Also rescale to preserve magnitude - shear[3] *= scalefac; - shear[4] *= scalefac; - shear[5] *= scalefac; - } - shear[3] += vrl1*dt; - shear[4] += vrl2*dt; - shear[5] += vrl3*dt; - } + // tangential forces = history + tangential velocity damping + fs1 = -k_tangential*history[thist0] - damp_tangential*vtr1; + fs2 = -k_tangential*history[thist1] - damp_tangential*vtr2; + fs3 = -k_tangential*history[thist2] - damp_tangential*vtr3; - if (rollingdamp == BRILLROLL) etaR = muR*fabs(Fne); - fr1 = -kR*shear[3] - etaR*vrl1; - fr2 = -kR*shear[4] - etaR*vrl2; - fr3 = -kR*shear[5] - etaR*vrl3; + // rescale frictional displacements and forces if needed + Fscrit = tangential_coeffs[2] * Fncrit; + fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); + if (fs > Fscrit) { + if (shrmag != 0.0) { + history[thist0] = -1.0/k_tangential*(Fscrit*fs1/fs + 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); + 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; + if (vrel != 0.0) Ft = MIN(Fne,fs) / vrel; + else Ft = 0.0; + fs1 = -Ft*vtr1; + fs2 = -Ft*vtr2; + fs3 = -Ft*vtr3; + } - // rescale frictional displacements and forces if needed - Frcrit = muR * fabs(Fne); - rollmag = sqrt(shear[3]*shear[3] + shear[4]*shear[4] + shear[5]*shear[5]); + //**************************************** + // Rolling resistance + //**************************************** - fr = sqrt(fr1*fr1 + fr2*fr2 + fr3*fr3); - if (fr > Frcrit) { - if (rollmag != 0.0) { - shear[3] = -1.0/kR*(Frcrit*fr1/fr + etaR*vrl1); - shear[4] = -1.0/kR*(Frcrit*fr2/fr + etaR*vrl2); - shear[5] = -1.0/kR*(Frcrit*fr3/fr + etaR*vrl3); - fr1 *= Frcrit/fr; - fr2 *= Frcrit/fr; - fr3 *= Frcrit/fr; - } else fr1 = fr2 = fr3 = 0.0; - } + if (roll_model != ROLL_NONE){ + relrot1 = omega[0]; + relrot2 = omega[1]; + relrot3 = omega[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) + vrl1 = Reff*(relrot2*nz - relrot3*ny); //- 0.5*((radj-radi)/radsum)*vtr1; + vrl2 = Reff*(relrot3*nx - relrot1*nz); //- 0.5*((radj-radi)/radsum)*vtr2; + vrl3 = Reff*(relrot1*ny - relrot2*nx); //- 0.5*((radj-radi)/radsum)*vtr3; + vrlmag = sqrt(vrl1*vrl1+vrl2*vrl2+vrl3*vrl3); + if (vrlmag != 0.0) vrlmaginv = 1.0/vrlmag; + else vrlmaginv = 0.0; - //**************************************** - // Twisting torque, including shear history effects - //**************************************** - magtwist = relrot1*nx + relrot2*ny + relrot3*nz; //Omega_T (eq 29 of Marshall) - shear[6] += magtwist*dt; - k_Q = 0.5*kt*a*a;; //eq 32 - eta_Q = 0.5*eta_T*a*a; - magtortwist = -k_Q*shear[6] - eta_Q*magtwist;//M_t torque (eq 30) + int rhist0 = roll_history_index; + int rhist1 = rhist0 + 1; + int rhist2 = rhist1 + 1; - signtwist = (magtwist > 0) - (magtwist < 0); - Mtcrit=TWOTHIRDS*a*Fscrit;//critical torque (eq 44) - if (fabs(magtortwist) > Mtcrit){ - shear[6] = 1.0/k_Q*(Mtcrit*signtwist - eta_Q*magtwist); - magtortwist = -Mtcrit * signtwist; //eq 34 - } + // Rolling displacement + rollmag = sqrt(history[rhist0]*history[rhist0] + + history[rhist1]*history[rhist1] + + history[rhist2]*history[rhist2]); - // Apply forces & torques + rolldotn = history[rhist0]*nx + history[rhist1]*ny + history[rhist2]*nz; - fx = nx*Fntot + fs1; - fy = ny*Fntot + fs2; - fz = nz*Fntot + fs3; + if (history_update){ + if (fabs(rolldotn) < EPSILON) rolldotn = 0; + if (rolldotn > 0){ //Rotate into tangential plane + scalefac = rollmag/(rollmag - rolldotn); + history[rhist0] -= rolldotn*nx; + history[rhist1] -= rolldotn*ny; + history[rhist2] -= rolldotn*nz; + //Also rescale to preserve magnitude + history[rhist0] *= scalefac; + history[rhist1] *= scalefac; + history[rhist2] *= scalefac; + } + history[rhist0] += vrl1*dt; + history[rhist1] += vrl2*dt; + history[rhist2] += vrl3*dt; + } - f[0] += fx; - f[1] += fy; - f[2] += fz; + 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; - tor1 = ny*fs3 - nz*fs2; - tor2 = nz*fs1 - nx*fs3; - tor3 = nx*fs2 - ny*fs1; + // rescale frictional displacements and forces if needed + Frcrit = roll_coeffs[2] * Fncrit; - torque[0] -= radi*tor1; - torque[1] -= radi*tor2; - torque[2] -= radi*tor3; + fr = sqrt(fr1*fr1 + fr2*fr2 + fr3*fr3); + if (fr > Frcrit) { + 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); + history[rhist2] = -1.0/k_roll*(Frcrit*fr3/fr + damp_roll*vrl3); + fr1 *= Frcrit/fr; + fr2 *= Frcrit/fr; + fr3 *= Frcrit/fr; + } else fr1 = fr2 = fr3 = 0.0; + } + } - tortwist1 = magtortwist * nx; - tortwist2 = magtortwist * ny; - tortwist3 = magtortwist * nz; + //**************************************** + // Twisting torque, including history effects + //**************************************** + if (twist_model != TWIST_NONE){ + magtwist = relrot1*nx + relrot2*ny + relrot3*nz; //Omega_T (eq 29 of Marshall) + if (twist_model == TWIST_MARSHALL){ + k_twist = 0.5*k_tangential*a*a;; //eq 32 of Marshall paper + damp_twist = 0.5*damp_tangential*a*a; + mu_twist = TWOTHIRDS*a*tangential_coeffs[2]; + } + else{ + k_twist = twist_coeffs[0]; + damp_twist = twist_coeffs[1]; + mu_twist = twist_coeffs[2]; + } + if (history_update){ + history[twist_history_index] += magtwist*dt; + } + magtortwist = -k_twist*history[twist_history_index] - damp_twist*magtwist;//M_t torque (eq 30) + signtwist = (magtwist > 0) - (magtwist < 0); + Mtcrit = mu_twist*Fncrit;//critical torque (eq 44) + if (fabs(magtortwist) > Mtcrit) { + history[twist_history_index] = 1.0/k_twist*(Mtcrit*signtwist - damp_twist*magtwist); + magtortwist = -Mtcrit * signtwist; //eq 34 + } + } + // Apply forces & torques - torque[0] += tortwist1; - torque[1] += tortwist2; - torque[2] += tortwist3; + fx = nx*Fntot + fs1; + fy = ny*Fntot + fs2; + fz = nz*Fntot + fs3; - torroll1 = R*(ny*fr3 - nz*fr2); //n cross fr - torroll2 = R*(nz*fr1 - nx*fr3); - torroll3 = R*(nx*fr2 - ny*fr1); + if (peratom_flag){ + contact[1] = fx; + contact[2] = fy; + contact[3] = fz; + } - torque[0] += torroll1; - torque[1] += torroll2; - torque[2] += torroll3; + f[0] += fx; + f[1] += fy; + f[2] += fz; + tor1 = ny*fs3 - nz*fs2; + tor2 = nz*fs1 - nx*fs3; + tor3 = nx*fs2 - ny*fs1; + + torque[0] -= radius*tor1; + torque[1] -= radius*tor2; + torque[2] -= radius*tor3; + + if (twist_model != TWIST_NONE){ + tortwist1 = magtortwist * nx; + tortwist2 = magtortwist * ny; + tortwist3 = magtortwist * nz; + + torque[0] += tortwist1; + torque[1] += tortwist2; + torque[2] += tortwist3; + } + + if (roll_model != ROLL_NONE){ + torroll1 = Reff*(ny*fr3 - nz*fr2); //n cross fr + torroll2 = Reff*(nz*fr1 - nx*fr3); + torroll3 = Reff*(nx*fr2 - ny*fr1); + + torque[0] += torroll1; + torque[1] += torroll2; + torque[2] += torroll3; + } } + /* ---------------------------------------------------------------------- memory usage of local atom-based arrays ------------------------------------------------------------------------- */ @@ -1124,7 +1383,7 @@ double FixWallGran::memory_usage() { int nmax = atom->nmax; double bytes = 0.0; - if (history) bytes += nmax*sheardim * sizeof(double); // shear history + if (use_history) bytes += nmax*size_history * sizeof(double); // shear history if (fix_rigid) bytes += nmax * sizeof(int); // mass_rigid if (peratom_flag) bytes += nmax*size_peratom_cols*sizeof(double); //store contacts return bytes; @@ -1136,7 +1395,7 @@ double FixWallGran::memory_usage() void FixWallGran::grow_arrays(int nmax) { - if (history) memory->grow(shearone,nmax,sheardim,"fix_wall_gran:shearone"); + if (use_history) memory->grow(history_one,nmax,size_history,"fix_wall_gran:history_one"); if (peratom_flag){ memory->grow(array_atom,nmax,size_peratom_cols,"fix_wall_gran:array_atom"); } @@ -1148,9 +1407,9 @@ void FixWallGran::grow_arrays(int nmax) void FixWallGran::copy_arrays(int i, int j, int /*delflag*/) { - if (history) - for (int m = 0; m < sheardim; m++) - shearone[j][m] = shearone[i][m]; + if (use_history) + for (int m = 0; m < size_history; m++) + history_one[j][m] = history_one[i][m]; if (peratom_flag){ for (int m = 0; m < size_peratom_cols; m++) array_atom[j][m] = array_atom[i][m]; @@ -1163,9 +1422,9 @@ void FixWallGran::copy_arrays(int i, int j, int /*delflag*/) void FixWallGran::set_arrays(int i) { - if (history) - for (int m = 0; m < sheardim; m++) - shearone[i][m] = 0; + if (use_history) + for (int m = 0; m < size_history; m++) + history_one[i][m] = 0; if (peratom_flag){ for (int m = 0; m < size_peratom_cols; m++) array_atom[i][m] = 0; @@ -1179,9 +1438,9 @@ void FixWallGran::set_arrays(int i) int FixWallGran::pack_exchange(int i, double *buf) { int n = 0; - if (history){ - for (int m = 0; m < sheardim; m++) - buf[n++] = shearone[i][m]; + if (use_history){ + for (int m = 0; m < size_history; m++) + buf[n++] = history_one[i][m]; } if (peratom_flag){ for (int m = 0; m < size_peratom_cols; m++) @@ -1197,9 +1456,9 @@ int FixWallGran::pack_exchange(int i, double *buf) int FixWallGran::unpack_exchange(int nlocal, double *buf) { int n = 0; - if (history){ - for (int m = 0; m < sheardim; m++) - shearone[nlocal][m] = buf[n++]; + if (use_history){ + for (int m = 0; m < size_history; m++) + history_one[nlocal][m] = buf[n++]; } if (peratom_flag){ for (int m = 0; m < size_peratom_cols; m++) @@ -1214,12 +1473,12 @@ int FixWallGran::unpack_exchange(int nlocal, double *buf) int FixWallGran::pack_restart(int i, double *buf) { - if (!history) return 0; + if (!use_history) return 0; int n = 0; - buf[n++] = sheardim + 1; - for (int m = 0; m < sheardim; m++) - buf[n++] = shearone[i][m]; + buf[n++] = size_history + 1; + for (int m = 0; m < size_history; m++) + buf[n++] = history_one[i][m]; return n; } @@ -1229,7 +1488,7 @@ int FixWallGran::pack_restart(int i, double *buf) void FixWallGran::unpack_restart(int nlocal, int nth) { - if (!history) return; + if (!use_history) return; double **extra = atom->extra; @@ -1239,8 +1498,8 @@ void FixWallGran::unpack_restart(int nlocal, int nth) for (int i = 0; i < nth; i++) m += static_cast (extra[nlocal][m]); m++; - for (int i = 0; i < sheardim; i++) - shearone[nlocal][i] = extra[nlocal][m++]; + for (int i = 0; i < size_history; i++) + history_one[nlocal][i] = extra[nlocal][m++]; } /* ---------------------------------------------------------------------- @@ -1249,8 +1508,8 @@ void FixWallGran::unpack_restart(int nlocal, int nth) int FixWallGran::maxsize_restart() { - if (!history) return 0; - return 1 + sheardim; + if (!use_history) return 0; + return 1 + size_history; } /* ---------------------------------------------------------------------- @@ -1259,8 +1518,8 @@ int FixWallGran::maxsize_restart() int FixWallGran::size_restart(int /*nlocal*/) { - if (!history) return 0; - return 1 + sheardim; + if (!use_history) return 0; + return 1 + size_history; } /* ---------------------------------------------------------------------- */ @@ -1270,3 +1529,12 @@ void FixWallGran::reset_dt() dt = update->dt; } +double FixWallGran::pulloff_distance(double radius){ + double coh, E, a, dist; + coh = normal_coeffs[3]; + E = normal_coeffs[0]*THREEQUARTERS; + a = cbrt(9*M_PI*coh*radius/(4*E)); + dist = a*a/radius - 2*sqrt(M_PI*coh*a/E); + return dist; +} + diff --git a/src/GRANULAR/fix_wall_gran.h b/src/GRANULAR/fix_wall_gran.h index 4212b96544..07c6c131cf 100644 --- a/src/GRANULAR/fix_wall_gran.h +++ b/src/GRANULAR/fix_wall_gran.h @@ -54,12 +54,11 @@ class FixWallGran : public Fix { void hertz_history(double, double, double, double, double *, double, double *, double *, double *, double *, double, double, double *, double *); - void dmt_rolling(double, double, double, double, double *, double, - double *, double *, double *, double *, double, double, - double *, double *); - // void jkr_rolling(double, double, double, double, double *, double, - // double *, double *, double *, double *, double, double, - // double *, double *); + void granular(double, double, double, double, double *, double, + double *, double *, double *, double *, double, double, + double *, double *); + + double pulloff_distance(double); protected: int wallstyle,wiggle,wshear,axis; @@ -67,23 +66,43 @@ class FixWallGran : public Fix { bigint time_origin; double kn,kt,gamman,gammat,xmu; - //For DMT/ROLLING - int normaldamp, rollingdamp; - double Emod, Gmod, alpha, Ecoh, kR, muR, etaR; + //For granular + //Model choices + int normal_model, damping_model; + int tangential_model, roll_model, twist_model; + int beyond_contact; + + //History flags + int normal_history, tangential_history, roll_history, twist_history; + + //Indices of history entries + int normal_history_index; + int tangential_history_index; + int roll_history_index; + int twist_history_index; + + //Material coefficients + double Emod, poiss, Gmod; + + //Contact model coefficients + double normal_coeffs[4]; + double tangential_coeffs[3]; + double roll_coeffs[3]; + double twist_coeffs[3]; double lo,hi,cylradius; double amplitude,period,omega,vshear; double dt; char *idregion; - int history; // if particle/wall interaction stores history - int shearupdate; // flag for whether shear history is updated - int sheardim; // # of shear history values per contact + int use_history; // if particle/wall interaction stores history + int history_update; // flag for whether shear history is updated + int size_history; // # of shear history values per contact // shear history for single contact per particle - double **shearone; + double **history_one; // rigid body masses for use in granular interactions diff --git a/src/GRANULAR/fix_wall_gran_region.cpp b/src/GRANULAR/fix_wall_gran_region.cpp index 17e16cb16b..95b34e0929 100644 --- a/src/GRANULAR/fix_wall_gran_region.cpp +++ b/src/GRANULAR/fix_wall_gran_region.cpp @@ -39,7 +39,8 @@ using namespace MathConst; // same as FixWallGran -enum{HOOKE,HOOKE_HISTORY,HERTZ_HISTORY,JKR_ROLLING,DMT_ROLLING}; +enum{HOOKE,HOOKE_HISTORY,HERTZ_HISTORY,GRANULAR}; +enum {NORMAL_HOOKE, NORMAL_HERTZ, HERTZ_MATERIAL, DMT, JKR}; #define BIG 1.0e20 @@ -47,7 +48,7 @@ enum{HOOKE,HOOKE_HISTORY,HERTZ_HISTORY,JKR_ROLLING,DMT_ROLLING}; FixWallGranRegion::FixWallGranRegion(LAMMPS *lmp, int narg, char **arg) : FixWallGran(lmp, narg, arg), region(NULL), region_style(NULL), ncontact(NULL), - walls(NULL), shearmany(NULL), c2r(NULL) + walls(NULL), history_many(NULL), c2r(NULL) { restart_global = 1; motion_resetflag = 0; @@ -66,17 +67,17 @@ FixWallGranRegion::FixWallGranRegion(LAMMPS *lmp, int narg, char **arg) : // re-allocate atom-based arrays with nshear // do not register with Atom class, since parent class did that - memory->destroy(shearone); - shearone = NULL; + memory->destroy(history_one); + history_one = NULL; ncontact = NULL; walls = NULL; - shearmany = NULL; + history_many = NULL; grow_arrays(atom->nmax); // initialize shear history as if particle is not touching region - if (history) { + if (use_history) { int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) ncontact[i] = 0; @@ -92,7 +93,7 @@ FixWallGranRegion::~FixWallGranRegion() memory->destroy(ncontact); memory->destroy(walls); - memory->destroy(shearmany); + memory->destroy(history_many); } /* ---------------------------------------------------------------------- */ @@ -138,8 +139,8 @@ void FixWallGranRegion::post_force(int /*vflag*/) // do not update shear history during setup - shearupdate = 1; - if (update->setupflag) shearupdate = 0; + history_update = 1; + if (update->setupflag) history_update = 0; // if just reneighbored: // update rigid body masses for owned atoms if using FixRigid @@ -188,7 +189,12 @@ void FixWallGranRegion::post_force(int /*vflag*/) if (mask[i] & groupbit) { if (!region->match(x[i][0],x[i][1],x[i][2])) continue; - nc = region->surface(x[i][0],x[i][1],x[i][2],radius[i]); + if (pairstyle == GRANULAR && normal_model == JKR){ + nc = region->surface(x[i][0],x[i][1],x[i][2],radius[i]+pulloff_distance(radius[i])); + } + else{ + nc = region->surface(x[i][0],x[i][1],x[i][2],radius[i]); + } if (nc > tmax) error->one(FLERR,"Too many wall/gran/region contacts for one particle"); @@ -198,7 +204,7 @@ void FixWallGranRegion::post_force(int /*vflag*/) // also set c2r[] = indices into region->contact[] for each of N contacts // process zero or one contact here, otherwise invoke update_contacts() - if (history) { + if (use_history) { if (nc == 0) { ncontact[i] = 0; continue; @@ -209,8 +215,8 @@ void FixWallGranRegion::post_force(int /*vflag*/) if (ncontact[i] == 0) { ncontact[i] = 1; walls[i][0] = iwall; - for (m = 0; m < sheardim; m++) - shearmany[i][0][m] = 0.0; + for (m = 0; m < size_history; m++) + history_many[i][0][m] = 0.0; } else if (ncontact[i] > 1 || iwall != walls[i][0]) update_contacts(i,nc); } else update_contacts(i,nc); @@ -224,13 +230,20 @@ void FixWallGranRegion::post_force(int /*vflag*/) rsq = region->contact[ic].r*region->contact[ic].r; + if (pairstyle == GRANULAR && normal_model == JKR){ + if (history_many[i][c2r[ic]][0] == 0.0 && rsq > radius[i]*radius[i]){ + for (m = 0; m < size_history; m++) + history_many[i][0][m] = 0.0; + continue; + } + } + dx = region->contact[ic].delx; dy = region->contact[ic].dely; dz = region->contact[ic].delz; if (regiondynamic) region->velocity_contact(vwall, x[i], ic); - // meff = effective mass of sphere // if I is part of rigid body, use body mass @@ -259,13 +272,13 @@ void FixWallGranRegion::post_force(int /*vflag*/) else if (pairstyle == HOOKE_HISTORY) hooke_history(rsq,dx,dy,dz,vwall,v[i],f[i], omega[i],torque[i],radius[i],meff, - shearmany[i][c2r[ic]], contact); + history_many[i][c2r[ic]], contact); else if (pairstyle == HERTZ_HISTORY) hertz_history(rsq,dx,dy,dz,vwall,region->contact[ic].radius, v[i],f[i],omega[i],torque[i], - radius[i],meff,shearmany[i][c2r[ic]], contact); - else if (pairstyle == DMT_ROLLING) - dmt_rolling(rsq,dx,dy,dz,vwall,region->contact[ic].radius, v[i],f[i],omega[i],torque[i], radius[i],meff,shearmany[i][c2r[ic]], contact); + radius[i],meff,history_many[i][c2r[ic]], contact); + else if (pairstyle == GRANULAR) + granular(rsq,dx,dy,dz,vwall,region->contact[ic].radius, v[i],f[i],omega[i],torque[i], radius[i],meff,history_many[i][c2r[ic]], contact); } } @@ -294,8 +307,8 @@ void FixWallGranRegion::update_contacts(int i, int nc) if (region->contact[m].iwall == walls[i][iold]) break; if (m >= nc) { ilast = ncontact[i]-1; - for (j = 0; j < sheardim; j++) - shearmany[i][iold][j] = shearmany[i][ilast][j]; + for (j = 0; j < size_history; j++) + history_many[i][iold][j] = history_many[i][ilast][j]; walls[i][iold] = walls[i][ilast]; ncontact[i]--; } else iold++; @@ -317,8 +330,8 @@ void FixWallGranRegion::update_contacts(int i, int nc) iadd = ncontact[i]; c2r[iadd] = inew; - for (j = 0; j < sheardim; j++) - shearmany[i][iadd][j] = 0.0; + for (j = 0; j < size_history; j++) + history_many[i][iadd][j] = 0.0; walls[i][iadd] = iwall; ncontact[i]++; } @@ -333,10 +346,10 @@ double FixWallGranRegion::memory_usage() { int nmax = atom->nmax; double bytes = 0.0; - if (history) { // shear history + if (use_history) { // shear history bytes += nmax * sizeof(int); // ncontact bytes += nmax*tmax * sizeof(int); // walls - bytes += nmax*tmax*sheardim * sizeof(double); // shearmany + bytes += nmax*tmax*size_history * sizeof(double); // history_many } if (fix_rigid) bytes += nmax * sizeof(int); // mass_rigid return bytes; @@ -348,10 +361,10 @@ double FixWallGranRegion::memory_usage() void FixWallGranRegion::grow_arrays(int nmax) { - if (history) { + if (use_history) { memory->grow(ncontact,nmax,"fix_wall_gran:ncontact"); memory->grow(walls,nmax,tmax,"fix_wall_gran:walls"); - memory->grow(shearmany,nmax,tmax,sheardim,"fix_wall_gran:shearmany"); + memory->grow(history_many,nmax,tmax,size_history,"fix_wall_gran:history_many"); } if (peratom_flag){ memory->grow(array_atom,nmax,size_peratom_cols,"fix_wall_gran:array_atom"); @@ -366,12 +379,12 @@ void FixWallGranRegion::copy_arrays(int i, int j, int /*delflag*/) { int m,n,iwall; - if (history){ + if (use_history){ n = ncontact[i]; for (iwall = 0; iwall < n; iwall++) { walls[j][iwall] = walls[i][iwall]; - for (m = 0; m < sheardim; m++) - shearmany[j][iwall][m] = shearmany[i][iwall][m]; + for (m = 0; m < size_history; m++) + history_many[j][iwall][m] = history_many[i][iwall][m]; } ncontact[j] = ncontact[i]; } @@ -388,7 +401,7 @@ void FixWallGranRegion::copy_arrays(int i, int j, int /*delflag*/) void FixWallGranRegion::set_arrays(int i) { - if (history) + if (use_history) ncontact[i] = 0; if (peratom_flag){ for (int m = 0; m < size_peratom_cols; m++) @@ -405,13 +418,13 @@ int FixWallGranRegion::pack_exchange(int i, double *buf) int m; int n = 0; - if (history){ + if (use_history){ int count = ncontact[i]; buf[n++] = ubuf(count).d; for (int iwall = 0; iwall < count; iwall++) { buf[n++] = ubuf(walls[i][iwall]).d; - for (m = 0; m < sheardim; m++) - buf[n++] = shearmany[i][iwall][m]; + for (m = 0; m < size_history; m++) + buf[n++] = history_many[i][iwall][m]; } } if (peratom_flag){ @@ -432,12 +445,12 @@ int FixWallGranRegion::unpack_exchange(int nlocal, double *buf) int n = 0; - if (history){ + if (use_history){ int count = ncontact[nlocal] = (int) ubuf(buf[n++]).i; for (int iwall = 0; iwall < count; iwall++) { walls[nlocal][iwall] = (int) ubuf(buf[n++]).i; - for (m = 0; m < sheardim; m++) - shearmany[nlocal][iwall][m] = buf[n++]; + for (m = 0; m < size_history; m++) + history_many[nlocal][iwall][m] = buf[n++]; } } if (peratom_flag){ @@ -456,7 +469,7 @@ int FixWallGranRegion::pack_restart(int i, double *buf) { int m; - if (!history) return 0; + if (!use_history) return 0; int n = 1; int count = ncontact[i]; @@ -464,8 +477,8 @@ int FixWallGranRegion::pack_restart(int i, double *buf) buf[n++] = ubuf(count).d; for (int iwall = 0; iwall < count; iwall++) { buf[n++] = ubuf(walls[i][iwall]).d; - for (m = 0; m < sheardim; m++) - buf[n++] = shearmany[i][iwall][m]; + for (m = 0; m < size_history; m++) + buf[n++] = history_many[i][iwall][m]; } buf[0] = n; return n; @@ -479,7 +492,7 @@ void FixWallGranRegion::unpack_restart(int nlocal, int nth) { int k; - if (!history) return; + if (!use_history) return; double **extra = atom->extra; @@ -492,8 +505,8 @@ void FixWallGranRegion::unpack_restart(int nlocal, int nth) int count = ncontact[nlocal] = (int) ubuf(extra[nlocal][m++]).i; for (int iwall = 0; iwall < count; iwall++) { walls[nlocal][iwall] = (int) ubuf(extra[nlocal][m++]).i; - for (k = 0; k < sheardim; k++) - shearmany[nlocal][iwall][k] = extra[nlocal][m++]; + for (k = 0; k < size_history; k++) + history_many[nlocal][iwall][k] = extra[nlocal][m++]; } } @@ -503,8 +516,8 @@ void FixWallGranRegion::unpack_restart(int nlocal, int nth) int FixWallGranRegion::maxsize_restart() { - if (!history) return 0; - return 2 + tmax*(sheardim+1); + if (!use_history) return 0; + return 2 + tmax*(size_history+1); } /* ---------------------------------------------------------------------- @@ -513,8 +526,8 @@ int FixWallGranRegion::maxsize_restart() int FixWallGranRegion::size_restart(int nlocal) { - if (!history) return 0; - return 2 + ncontact[nlocal]*(sheardim+1); + if (!use_history) return 0; + return 2 + ncontact[nlocal]*(size_history+1); } /* ---------------------------------------------------------------------- diff --git a/src/GRANULAR/fix_wall_gran_region.h b/src/GRANULAR/fix_wall_gran_region.h index 8d1b6d533a..fd40e27e4c 100644 --- a/src/GRANULAR/fix_wall_gran_region.h +++ b/src/GRANULAR/fix_wall_gran_region.h @@ -54,7 +54,7 @@ class FixWallGranRegion : public FixWallGran { int tmax; // max # of region walls one particle can touch int *ncontact; // # of shear contacts per particle int **walls; // which wall each contact is with - double ***shearmany; // shear history per particle per contact + double ***history_many; // history per particle per contact int *c2r; // contact to region mapping // c2r[i] = index of Ith contact in // region-contact[] list of contacts diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index 5631240fea..ac0b668854 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -45,6 +45,7 @@ using namespace MathConst; #define SIXROOT6 14.69693845669906728801 // 6*sqrt(6) #define INVROOT6 0.40824829046386307274 // 1/sqrt(6) #define FOURTHIRDS 1.333333333333333 // 4/3 +#define THREEQUARTERS 0.75 // 3/4 #define TWOPI 6.28318530717959 // 2*PI #define EPSILON 1e-10 @@ -63,7 +64,7 @@ PairGranular::PairGranular(LAMMPS *lmp) : Pair(lmp) no_virial_fdotr_compute = 1; fix_history = NULL; - single_extra = 9; + single_extra = 12; svector = new double[single_extra]; neighprev = 0; @@ -238,10 +239,11 @@ void PairGranular::compute(int eflag, int vflag) radsum = radi + radj; E = normal_coeffs[itype][jtype][0]; - Reff = radi*radj/(radi+radj); + Reff = radi*radj/radsum; touchflag = false; if (normal_model[itype][jtype] == JKR){ + E *= THREEQUARTERS; if (touch[jj]){ R2 = Reff*Reff; coh = normal_coeffs[itype][jtype][3]; @@ -319,17 +321,17 @@ void PairGranular::compute(int eflag, int vflag) sqrt3 = MAX(0, 4*dR - t5 + SIXROOT6*coh*M_PI*R2/(E*t6)); a = INVROOT6*(t6 + sqrt(sqrt3)); a2 = a*a; - knfac = FOURTHIRDS*E*a; + knfac = normal_coeffs[itype][jtype][0]*a; Fne = knfac*a2/Reff - TWOPI*a2*sqrt(4*coh*E/(M_PI*a)); } else{ knfac = E; //Hooke + Fne = knfac*delta; a = sqrt(dR); if (normal_model[itype][jtype] != HOOKE){ Fne *= a; knfac *= a; } - Fne = knfac*delta; if (normal_model[itype][jtype] == DMT) Fne -= 4*MY_PI*normal_coeffs[itype][jtype][3]*Reff; } @@ -711,9 +713,9 @@ void PairGranular::coeff(int narg, char **arg) } else if (strcmp(arg[iarg], "hertz/material") == 0){ int num_coeffs = 3; - if (iarg + num_coeffs >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for Hertz option"); - normal_model_one = HERTZ; - normal_coeffs_one[0] = force->numeric(FLERR,arg[iarg+1])*FOURTHIRDS; //E + if (iarg + num_coeffs >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for Hertz/material option"); + normal_model_one = HERTZ_MATERIAL; + normal_coeffs_one[0] = force->numeric(FLERR,arg[iarg+1]); //E normal_coeffs_one[1] = force->numeric(FLERR,arg[iarg+2]); //damping normal_coeffs_one[2] = force->numeric(FLERR,arg[iarg+3]); //Poisson's ratio iarg += num_coeffs+1; @@ -721,10 +723,10 @@ void PairGranular::coeff(int narg, char **arg) else if (strcmp(arg[iarg], "dmt") == 0){ if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for Hertz option"); normal_model_one = DMT; - normal_coeffs_one[0] = force->numeric(FLERR,arg[iarg+1])*FOURTHIRDS; //E + normal_coeffs_one[0] = force->numeric(FLERR,arg[iarg+1]); //E normal_coeffs_one[1] = force->numeric(FLERR,arg[iarg+2]); //damping normal_coeffs_one[2] = force->numeric(FLERR,arg[iarg+3]); //Poisson's ratio - normal_coeffs_one[3] = force->numeric(FLERR,arg[iarg+3]); //cohesion + normal_coeffs_one[3] = force->numeric(FLERR,arg[iarg+4]); //cohesion iarg += 5; } else if (strcmp(arg[iarg], "jkr") == 0){ @@ -755,21 +757,27 @@ void PairGranular::coeff(int narg, char **arg) iarg += 1; } else if (strcmp(arg[iarg], "tangential") == 0){ - if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for tangential model"); + if (iarg + 1 >= narg) error->all(FLERR,"Illegal pair_coeff command, must specify tangential model after 'tangential' keyword"); if (strcmp(arg[iarg+1], "linear_nohistory") == 0){ + if (iarg + 3 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for tangential model"); tangential_model_one = TANGENTIAL_NOHISTORY; + tangential_coeffs_one[0] = 0; + tangential_coeffs_one[1] = force->numeric(FLERR,arg[iarg+2]); //gammat + tangential_coeffs_one[2] = force->numeric(FLERR,arg[iarg+3]); //friction coeff. + iarg += 4; } else if (strcmp(arg[iarg+1], "linear_history") == 0){ + if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for tangential model"); tangential_model_one = TANGENTIAL_HISTORY; 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 + tangential_coeffs_one[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff. + iarg += 5; } else{ error->all(FLERR, "Illegal pair_coeff command, tangential model not recognized"); } - 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; } else if (strcmp(arg[iarg], "rolling") == 0){ if (iarg + 1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters"); @@ -842,7 +850,7 @@ void PairGranular::coeff(int narg, char **arg) if (normal_model_one != HERTZ && normal_model_one != HOOKE){ Emod[i][j] = Emod[j][i] = normal_coeffs_one[0]; poiss[i][j] = poiss[j][i] = normal_coeffs_one[2]; - normal_coeffs[i][j][0] = normal_coeffs[j][i][0] = mix_stiffnessE(Emod[i][j], Emod[i][j], poiss[i][j], poiss[i][j]); + normal_coeffs[i][j][0] = normal_coeffs[j][i][0] = FOURTHIRDS*mix_stiffnessE(Emod[i][j], Emod[i][j], poiss[i][j], poiss[i][j]); } else{ normal_coeffs[i][j][0] = normal_coeffs[j][i][0] = normal_coeffs_one[0]; @@ -1068,14 +1076,13 @@ double PairGranular::init_one(int i, int j) if (((maxrad_dynamic[i] > 0.0) && (maxrad_dynamic[j] > 0.0)) || ((maxrad_dynamic[i] > 0.0) && (maxrad_frozen[j] > 0.0)) || ((maxrad_frozen[i] > 0.0) && (maxrad_dynamic[j] > 0.0))) { // radius info about both i and j exist + cutoff = maxrad_dynamic[i]+maxrad_dynamic[j]; + pulloff = 0.0; if (normal_model[i][j] == JKR){ pulloff = pulloff_distance(maxrad_dynamic[i], maxrad_dynamic[j], i, j); cutoff += pulloff; } - else{ - pulloff = 0; - } if (normal_model[i][j] == JKR) pulloff = pulloff_distance(maxrad_frozen[i], maxrad_dynamic[j], i, j); @@ -1224,10 +1231,12 @@ double PairGranular::single(int i, int j, int itype, int jtype, radi = radius[i]; radj = radius[j]; radsum = radi + radj; - Reff = radi*radj/(radi+radj); + Reff = radi*radj/radsum; bool touchflag; + E = normal_coeffs[itype][jtype][0]; if (normal_model[itype][jtype] == JKR){ + E *= THREEQUARTERS; R2 = Reff*Reff; coh = normal_coeffs[itype][jtype][3]; a = cbrt(9.0*M_PI*coh*R2/(4*E)); @@ -1333,7 +1342,7 @@ double PairGranular::single(int i, int j, int itype, int jtype, sqrt3 = MAX(0, 4*dR - t5 + SIXROOT6*coh*M_PI*R2/(E*t6)); a = INVROOT6*(t6 + sqrt(sqrt3)); a2 = a*a; - knfac = FOURTHIRDS*E*a; + knfac = normal_coeffs[itype][jtype][0]*a; Fne = knfac*a2/Reff - TWOPI*a2*sqrt(4*coh*E/(M_PI*a)); } else{ @@ -1348,7 +1357,6 @@ double PairGranular::single(int i, int j, int itype, int jtype, Fne -= 4*MY_PI*normal_coeffs[itype][jtype][3]*Reff; } - //Consider restricting Hooke to only have 'velocity' as an option for damping? if (damping_model[itype][jtype] == VELOCITY){ damp_normal = normal_coeffs[itype][jtype][1]; } @@ -1536,6 +1544,9 @@ double PairGranular::single(int i, int j, int itype, int jtype, svector[6] = fr3; svector[7] = fr; svector[8] = magtortwist; + svector[9] = delx; + svector[10] = dely; + svector[11] = delz; return 0.0; } @@ -1614,7 +1625,7 @@ double PairGranular::pulloff_distance(double radi, double radj, int itype, int j Reff = radi*radj/(radi+radj); if (Reff <= 0) return 0; coh = normal_coeffs[itype][itype][3]; - E = normal_coeffs[itype][jtype][0]; + E = normal_coeffs[itype][jtype][0]*THREEQUARTERS; a = cbrt(9*M_PI*coh*Reff/(4*E)); return a*a/Reff - 2*sqrt(M_PI*coh*a/E); } diff --git a/src/GRANULAR/pair_granular.h b/src/GRANULAR/pair_granular.h index 625ff17c72..7bce3831f1 100644 --- a/src/GRANULAR/pair_granular.h +++ b/src/GRANULAR/pair_granular.h @@ -65,9 +65,9 @@ public: private: int size_history; - //Models choices + //Model choices int **normal_model, **damping_model; - double **tangential_model, **roll_model, **twist_model; + int **tangential_model, **roll_model, **twist_model; //History flags int normal_history, tangential_history, roll_history, twist_history;