forked from lijiext/lammps
change source code format style to be more like other LAMMPS sources
This commit is contained in:
Axel Kohlmeyer
parent
05a5ecd4d4
commit
4cd0ea61f2
|
|
@ -147,18 +147,18 @@ void PairGranular::compute(int eflag, int vflag)
|
|||
double mi,mj,meff;
|
||||
double relrot1,relrot2,relrot3,vrl1,vrl2,vrl3;
|
||||
|
||||
//For JKR
|
||||
// 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;
|
||||
|
||||
//Rolling
|
||||
// Rolling
|
||||
double k_roll, damp_roll;
|
||||
double torroll1, torroll2, torroll3;
|
||||
double rollmag, rolldotn, scalefac;
|
||||
double fr, fr1, fr2, fr3;
|
||||
|
||||
//Twisting
|
||||
// Twisting
|
||||
double k_twist, damp_twist, mu_twist;
|
||||
double signtwist, magtwist, magtortwist, Mtcrit;
|
||||
double tortwist1, tortwist2, tortwist3;
|
||||
|
|
@ -180,7 +180,7 @@ void PairGranular::compute(int eflag, int vflag)
|
|||
// body[i] = which body atom I is in, -1 if none
|
||||
// mass_body = mass of each rigid body
|
||||
|
||||
if (fix_rigid && neighbor->ago == 0){
|
||||
if (fix_rigid && neighbor->ago == 0) {
|
||||
int tmp;
|
||||
int *body = (int *) fix_rigid->extract("body",tmp);
|
||||
double *mass_body = (double *) fix_rigid->extract("masstotal",tmp);
|
||||
|
|
@ -211,7 +211,7 @@ void PairGranular::compute(int eflag, int vflag)
|
|||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
if (use_history){
|
||||
if (use_history) {
|
||||
firsttouch = fix_history->firstflag;
|
||||
firsthistory = fix_history->firstvalue;
|
||||
}
|
||||
|
|
@ -224,14 +224,14 @@ void PairGranular::compute(int eflag, int vflag)
|
|||
ztmp = x[i][2];
|
||||
itype = type[i];
|
||||
radi = radius[i];
|
||||
if (use_history){
|
||||
if (use_history) {
|
||||
touch = firsttouch[i];
|
||||
allhistory = firsthistory[i];
|
||||
}
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
|
||||
for (jj = 0; jj < jnum; jj++){
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
|
|
@ -247,33 +247,30 @@ void PairGranular::compute(int eflag, int vflag)
|
|||
Reff = radi*radj/radsum;
|
||||
touchflag = false;
|
||||
|
||||
if (normal_model[itype][jtype] == JKR){
|
||||
if (normal_model[itype][jtype] == JKR) {
|
||||
E *= THREEQUARTERS;
|
||||
if (touch[jj]){
|
||||
if (touch[jj]) {
|
||||
R2 = Reff*Reff;
|
||||
coh = normal_coeffs[itype][jtype][3];
|
||||
a = cbrt(9.0*M_PI*coh*R2/(4*E));
|
||||
delta_pulloff = a*a/Reff - 2*sqrt(M_PI*coh*a/E);
|
||||
dist_pulloff = radsum-delta_pulloff;
|
||||
touchflag = (rsq < dist_pulloff*dist_pulloff);
|
||||
}
|
||||
else{
|
||||
} else {
|
||||
touchflag = (rsq < radsum*radsum);
|
||||
}
|
||||
}
|
||||
else{
|
||||
} else {
|
||||
touchflag = (rsq < radsum*radsum);
|
||||
}
|
||||
|
||||
if (!touchflag){
|
||||
if (!touchflag) {
|
||||
// unset non-touching neighbors
|
||||
if (use_history){
|
||||
if (use_history) {
|
||||
touch[jj] = 0;
|
||||
history = &allhistory[size_history*jj];
|
||||
for (int k = 0; k < size_history; k++) history[k] = 0.0;
|
||||
}
|
||||
}
|
||||
else{
|
||||
} else {
|
||||
r = sqrt(rsq);
|
||||
rinv = 1.0/r;
|
||||
|
||||
|
|
@ -311,7 +308,7 @@ void PairGranular::compute(int eflag, int vflag)
|
|||
|
||||
delta = radsum - r;
|
||||
dR = delta*Reff;
|
||||
if (normal_model[itype][jtype] == JKR){
|
||||
if (normal_model[itype][jtype] == JKR) {
|
||||
touch[jj] = 1;
|
||||
R2=Reff*Reff;
|
||||
coh = normal_coeffs[itype][jtype][3];
|
||||
|
|
@ -330,12 +327,11 @@ void PairGranular::compute(int eflag, int vflag)
|
|||
a2 = a*a;
|
||||
knfac = normal_coeffs[itype][jtype][0]*a;
|
||||
Fne = knfac*a2/Reff - TWOPI*a2*sqrt(4*coh*E/(M_PI*a));
|
||||
}
|
||||
else{
|
||||
} else {
|
||||
knfac = E; //Hooke
|
||||
Fne = knfac*delta;
|
||||
a = sqrt(dR);
|
||||
if (normal_model[itype][jtype] != HOOKE){
|
||||
if (normal_model[itype][jtype] != HOOKE) {
|
||||
Fne *= a;
|
||||
knfac *= a;
|
||||
}
|
||||
|
|
@ -344,13 +340,11 @@ void PairGranular::compute(int eflag, int vflag)
|
|||
}
|
||||
|
||||
//Consider restricting Hooke to only have 'velocity' as an option for damping?
|
||||
if (damping_model[itype][jtype] == VELOCITY){
|
||||
if (damping_model[itype][jtype] == VELOCITY) {
|
||||
damp_normal = 1;
|
||||
}
|
||||
else if (damping_model[itype][jtype] == VISCOELASTIC){
|
||||
} else if (damping_model[itype][jtype] == VISCOELASTIC) {
|
||||
damp_normal = a*meff;
|
||||
}
|
||||
else if (damping_model[itype][jtype] == TSUJI){
|
||||
} else if (damping_model[itype][jtype] == TSUJI) {
|
||||
damp_normal = sqrt(meff*knfac);
|
||||
}
|
||||
|
||||
|
|
@ -381,20 +375,18 @@ void PairGranular::compute(int eflag, int vflag)
|
|||
vrel = sqrt(vrel);
|
||||
|
||||
// If any history is needed:
|
||||
if (use_history){
|
||||
if (use_history) {
|
||||
touch[jj] = 1;
|
||||
history = &allhistory[size_history*jj];
|
||||
}
|
||||
|
||||
if (normal_model[itype][jtype] == JKR){
|
||||
if (normal_model[itype][jtype] == JKR) {
|
||||
F_pulloff = 3*M_PI*coh*Reff;
|
||||
Fncrit = fabs(Fne + 2*F_pulloff);
|
||||
}
|
||||
else if (normal_model[itype][jtype] == DMT){
|
||||
} else if (normal_model[itype][jtype] == DMT) {
|
||||
F_pulloff = 4*M_PI*coh*Reff;
|
||||
Fncrit = fabs(Fne + 2*F_pulloff);
|
||||
}
|
||||
else{
|
||||
} else {
|
||||
Fncrit = fabs(Fntot);
|
||||
}
|
||||
|
||||
|
|
@ -404,13 +396,12 @@ void PairGranular::compute(int eflag, int vflag)
|
|||
k_tangential = tangential_coeffs[itype][jtype][0];
|
||||
damp_tangential = tangential_coeffs[itype][jtype][1]*damp_normal_prefactor;
|
||||
|
||||
if (tangential_history){
|
||||
if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN){
|
||||
if (tangential_history) {
|
||||
if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN) {
|
||||
k_tangential *= a;
|
||||
}
|
||||
else if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE){
|
||||
} else if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) {
|
||||
k_tangential *= a;
|
||||
if (a < history[3]){ //On unloading, rescale the shear displacements
|
||||
if (a < history[3]) { //On unloading, rescale the shear displacements
|
||||
double factor = a/history[3];
|
||||
history[0] *= factor;
|
||||
history[1] *= factor;
|
||||
|
|
@ -422,7 +413,7 @@ void PairGranular::compute(int eflag, int vflag)
|
|||
if (historyupdate) {
|
||||
rsht = history[0]*nx + history[1]*ny + history[2]*nz;
|
||||
if (fabs(rsht) < EPSILON) rsht = 0;
|
||||
if (rsht > 0){
|
||||
if (rsht > 0) {
|
||||
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!
|
||||
|
|
@ -461,8 +452,7 @@ void PairGranular::compute(int eflag, int vflag)
|
|||
fs3 *= Fscrit/fs;
|
||||
} else fs1 = fs2 = fs3 = 0.0;
|
||||
}
|
||||
}
|
||||
else{ //Classic pair gran/hooke (no history)
|
||||
} 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;
|
||||
|
|
@ -475,7 +465,7 @@ void PairGranular::compute(int eflag, int vflag)
|
|||
// Rolling resistance
|
||||
//****************************************
|
||||
|
||||
if (roll_model[itype][jtype] != ROLL_NONE){
|
||||
if (roll_model[itype][jtype] != ROLL_NONE) {
|
||||
relrot1 = omega[i][0] - omega[j][0];
|
||||
relrot2 = omega[i][1] - omega[j][1];
|
||||
relrot3 = omega[i][2] - omega[j][2];
|
||||
|
|
@ -492,9 +482,9 @@ void PairGranular::compute(int eflag, int vflag)
|
|||
int rhist2 = rhist1 + 1;
|
||||
|
||||
rolldotn = history[rhist0]*nx + history[rhist1]*ny + history[rhist2]*nz;
|
||||
if (historyupdate){
|
||||
if (historyupdate) {
|
||||
if (fabs(rolldotn) < EPSILON) rolldotn = 0;
|
||||
if (rolldotn > 0){ //Rotate into tangential plane
|
||||
if (rolldotn > 0) { //Rotate into tangential plane
|
||||
rollmag = sqrt(history[rhist0]*history[rhist0] +
|
||||
history[rhist1]*history[rhist1] +
|
||||
history[rhist2]*history[rhist2]);
|
||||
|
|
@ -540,19 +530,18 @@ void PairGranular::compute(int eflag, int vflag)
|
|||
//****************************************
|
||||
// Twisting torque, including history effects
|
||||
//****************************************
|
||||
if (twist_model[itype][jtype] != TWIST_NONE){
|
||||
if (twist_model[itype][jtype] != TWIST_NONE) {
|
||||
magtwist = relrot1*nx + relrot2*ny + relrot3*nz; //Omega_T (eq 29 of Marshall)
|
||||
if (twist_model[itype][jtype] == TWIST_MARSHALL){
|
||||
if (twist_model[itype][jtype] == 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[itype][jtype][2];
|
||||
}
|
||||
else{
|
||||
} else {
|
||||
k_twist = twist_coeffs[itype][jtype][0];
|
||||
damp_twist = twist_coeffs[itype][jtype][1];
|
||||
mu_twist = twist_coeffs[itype][jtype][2];
|
||||
}
|
||||
if (historyupdate){
|
||||
if (historyupdate) {
|
||||
history[twist_history_index] += magtwist*dt;
|
||||
}
|
||||
magtortwist = -k_twist*history[twist_history_index] - damp_twist*magtwist;//M_t torque (eq 30)
|
||||
|
|
@ -582,7 +571,7 @@ void PairGranular::compute(int eflag, int vflag)
|
|||
torque[i][1] -= dist_to_contact*tor2;
|
||||
torque[i][2] -= dist_to_contact*tor3;
|
||||
|
||||
if (twist_model[itype][jtype] != TWIST_NONE){
|
||||
if (twist_model[itype][jtype] != TWIST_NONE) {
|
||||
tortwist1 = magtortwist * nx;
|
||||
tortwist2 = magtortwist * ny;
|
||||
tortwist3 = magtortwist * nz;
|
||||
|
|
@ -592,7 +581,7 @@ void PairGranular::compute(int eflag, int vflag)
|
|||
torque[i][2] += tortwist3;
|
||||
}
|
||||
|
||||
if (roll_model[itype][jtype] != ROLL_NONE){
|
||||
if (roll_model[itype][jtype] != ROLL_NONE) {
|
||||
torroll1 = Reff*(ny*fr3 - nz*fr2); //n cross fr
|
||||
torroll2 = Reff*(nz*fr1 - nx*fr3);
|
||||
torroll3 = Reff*(nx*fr2 - ny*fr1);
|
||||
|
|
@ -612,12 +601,12 @@ void PairGranular::compute(int eflag, int vflag)
|
|||
torque[j][1] -= dist_to_contact*tor2;
|
||||
torque[j][2] -= dist_to_contact*tor3;
|
||||
|
||||
if (twist_model[itype][jtype] != TWIST_NONE){
|
||||
if (twist_model[itype][jtype] != TWIST_NONE) {
|
||||
torque[j][0] -= tortwist1;
|
||||
torque[j][1] -= tortwist2;
|
||||
torque[j][2] -= tortwist3;
|
||||
}
|
||||
if (roll_model[itype][jtype] != ROLL_NONE){
|
||||
if (roll_model[itype][jtype] != ROLL_NONE) {
|
||||
torque[j][0] -= torroll1;
|
||||
torque[j][1] -= torroll2;
|
||||
torque[j][2] -= torroll3;
|
||||
|
|
@ -673,10 +662,9 @@ void PairGranular::allocate()
|
|||
|
||||
void PairGranular::settings(int narg, char **arg)
|
||||
{
|
||||
if (narg == 1){
|
||||
if (narg == 1) {
|
||||
cutoff_global = force->numeric(FLERR,arg[0]);
|
||||
}
|
||||
else{
|
||||
} else {
|
||||
cutoff_global = -1; //Will be set based on particle sizes, model choice
|
||||
}
|
||||
|
||||
|
|
@ -715,23 +703,21 @@ void PairGranular::coeff(int narg, char **arg)
|
|||
damping_model_one = VISCOELASTIC;
|
||||
|
||||
int iarg = 2;
|
||||
while (iarg < narg){
|
||||
if (strcmp(arg[iarg], "hooke") == 0){
|
||||
while (iarg < narg) {
|
||||
if (strcmp(arg[iarg], "hooke") == 0) {
|
||||
if (iarg + 2 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for Hooke option");
|
||||
normal_model_one = HOOKE;
|
||||
normal_coeffs_one[0] = force->numeric(FLERR,arg[iarg+1]); //kn
|
||||
normal_coeffs_one[1] = force->numeric(FLERR,arg[iarg+2]); //damping
|
||||
iarg += 3;
|
||||
}
|
||||
else if (strcmp(arg[iarg], "hertz") == 0){
|
||||
} else if (strcmp(arg[iarg], "hertz") == 0) {
|
||||
int num_coeffs = 2;
|
||||
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]); //kn
|
||||
normal_coeffs_one[1] = force->numeric(FLERR,arg[iarg+2]); //damping
|
||||
iarg += num_coeffs+1;
|
||||
}
|
||||
else if (strcmp(arg[iarg], "hertz/material") == 0){
|
||||
} 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/material option");
|
||||
normal_model_one = HERTZ_MATERIAL;
|
||||
|
|
@ -739,8 +725,7 @@ void PairGranular::coeff(int narg, char **arg)
|
|||
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;
|
||||
}
|
||||
else if (strcmp(arg[iarg], "dmt") == 0){
|
||||
} 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]); //E
|
||||
|
|
@ -748,8 +733,7 @@ void PairGranular::coeff(int narg, char **arg)
|
|||
normal_coeffs_one[2] = force->numeric(FLERR,arg[iarg+3]); //Poisson's ratio
|
||||
normal_coeffs_one[3] = force->numeric(FLERR,arg[iarg+4]); //cohesion
|
||||
iarg += 5;
|
||||
}
|
||||
else if (strcmp(arg[iarg], "jkr") == 0){
|
||||
} else if (strcmp(arg[iarg], "jkr") == 0) {
|
||||
if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for JKR option");
|
||||
beyond_contact = 1;
|
||||
normal_model_one = JKR;
|
||||
|
|
@ -758,67 +742,57 @@ void PairGranular::coeff(int narg, char **arg)
|
|||
normal_coeffs_one[2] = force->numeric(FLERR,arg[iarg+3]); //Poisson's ratio
|
||||
normal_coeffs_one[3] = force->numeric(FLERR,arg[iarg+4]); //cohesion
|
||||
iarg += 5;
|
||||
}
|
||||
else if (strcmp(arg[iarg], "damping") == 0){
|
||||
} else if (strcmp(arg[iarg], "damping") == 0) {
|
||||
if (iarg+1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters provided for damping model");
|
||||
if (strcmp(arg[iarg+1], "velocity") == 0){
|
||||
if (strcmp(arg[iarg+1], "velocity") == 0) {
|
||||
damping_model_one = VELOCITY;
|
||||
iarg += 1;
|
||||
}
|
||||
else if (strcmp(arg[iarg+1], "viscoelastic") == 0){
|
||||
} else if (strcmp(arg[iarg+1], "viscoelastic") == 0) {
|
||||
damping_model_one = VISCOELASTIC;
|
||||
iarg += 1;
|
||||
}
|
||||
else if (strcmp(arg[iarg+1], "tsuji") == 0){
|
||||
} else if (strcmp(arg[iarg+1], "tsuji") == 0) {
|
||||
damping_model_one = TSUJI;
|
||||
iarg += 1;
|
||||
}
|
||||
else error->all(FLERR, "Illegal pair_coeff command, unrecognized damping model");
|
||||
} else error->all(FLERR, "Illegal pair_coeff command, unrecognized damping model");
|
||||
iarg += 1;
|
||||
}
|
||||
else if (strcmp(arg[iarg], "tangential") == 0){
|
||||
} 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 (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) ||
|
||||
} else if ((strcmp(arg[iarg+1], "linear_history") == 0) ||
|
||||
(strcmp(arg[iarg+1], "mindlin") == 0) ||
|
||||
(strcmp(arg[iarg+1], "mindlin_rescale") == 0)){
|
||||
(strcmp(arg[iarg+1], "mindlin_rescale") == 0)) {
|
||||
if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for tangential model");
|
||||
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;
|
||||
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){
|
||||
(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{
|
||||
} 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;
|
||||
}
|
||||
else{
|
||||
} else {
|
||||
error->all(FLERR, "Illegal pair_coeff command, tangential model not recognized");
|
||||
}
|
||||
}
|
||||
else if (strcmp(arg[iarg], "rolling") == 0){
|
||||
} else if (strcmp(arg[iarg], "rolling") == 0) {
|
||||
if (iarg + 1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters");
|
||||
if (strcmp(arg[iarg+1], "none") == 0){
|
||||
if (strcmp(arg[iarg+1], "none") == 0) {
|
||||
roll_model_one = ROLL_NONE;
|
||||
iarg += 2;
|
||||
}
|
||||
else if (strcmp(arg[iarg+1], "sds") == 0){
|
||||
} else if (strcmp(arg[iarg+1], "sds") == 0) {
|
||||
if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for rolling model");
|
||||
roll_model_one = ROLL_SDS;
|
||||
roll_history = 1;
|
||||
|
|
@ -826,23 +800,19 @@ void PairGranular::coeff(int narg, char **arg)
|
|||
roll_coeffs_one[1] = force->numeric(FLERR,arg[iarg+3]); //gammaR
|
||||
roll_coeffs_one[2] = force->numeric(FLERR,arg[iarg+4]); //rolling friction coeff.
|
||||
iarg += 5;
|
||||
}
|
||||
else{
|
||||
} else {
|
||||
error->all(FLERR, "Illegal pair_coeff command, rolling friction model not recognized");
|
||||
}
|
||||
}
|
||||
else if (strcmp(arg[iarg], "twisting") == 0){
|
||||
} else if (strcmp(arg[iarg], "twisting") == 0) {
|
||||
if (iarg + 1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters");
|
||||
if (strcmp(arg[iarg+1], "none") == 0){
|
||||
if (strcmp(arg[iarg+1], "none") == 0) {
|
||||
twist_model_one = TWIST_NONE;
|
||||
iarg += 2;
|
||||
}
|
||||
else if (strcmp(arg[iarg+1], "marshall") == 0){
|
||||
} else if (strcmp(arg[iarg+1], "marshall") == 0) {
|
||||
twist_model_one = TWIST_MARSHALL;
|
||||
twist_history = 1;
|
||||
iarg += 2;
|
||||
}
|
||||
else if (strcmp(arg[iarg+1], "sds") == 0){
|
||||
} else if (strcmp(arg[iarg+1], "sds") == 0) {
|
||||
if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for twist model");
|
||||
twist_model_one = TWIST_SDS;
|
||||
twist_history = 1;
|
||||
|
|
@ -850,16 +820,13 @@ void PairGranular::coeff(int narg, char **arg)
|
|||
twist_coeffs_one[1] = force->numeric(FLERR,arg[iarg+3]); //gammat
|
||||
twist_coeffs_one[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff.
|
||||
iarg += 5;
|
||||
}
|
||||
else{
|
||||
} else {
|
||||
error->all(FLERR, "Illegal pair_coeff command, twisting friction model not recognized");
|
||||
}
|
||||
}
|
||||
else if (strcmp(arg[iarg], "cutoff") == 0){
|
||||
} else if (strcmp(arg[iarg], "cutoff") == 0) {
|
||||
if (iarg + 1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters");
|
||||
cutoff_one = force->numeric(FLERR,arg[iarg+1]);
|
||||
}
|
||||
else error->all(FLERR, "Illegal pair coeff command");
|
||||
} else error->all(FLERR, "Illegal pair coeff command");
|
||||
}
|
||||
|
||||
//It is an error not to specify normal or tangential model
|
||||
|
|
@ -867,25 +834,23 @@ void PairGranular::coeff(int narg, char **arg)
|
|||
|
||||
int count = 0;
|
||||
double damp;
|
||||
if (damping_model_one == TSUJI){
|
||||
if (damping_model_one == TSUJI) {
|
||||
double cor;
|
||||
cor = normal_coeffs_one[1];
|
||||
damp = 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);
|
||||
}
|
||||
else damp = normal_coeffs_one[1];
|
||||
} else damp = normal_coeffs_one[1];
|
||||
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
for (int j = MAX(jlo,i); j <= jhi; j++) {
|
||||
normal_model[i][j] = normal_model[j][i] = normal_model_one;
|
||||
normal_coeffs[i][j][1] = normal_coeffs[j][i][1] = damp;
|
||||
if (normal_model_one != HERTZ && normal_model_one != HOOKE){
|
||||
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] = FOURTHIRDS*mix_stiffnessE(Emod[i][j], Emod[i][j], poiss[i][j], poiss[i][j]);
|
||||
}
|
||||
else{
|
||||
} else {
|
||||
normal_coeffs[i][j][0] = normal_coeffs[j][i][0] = normal_coeffs_one[0];
|
||||
}
|
||||
if ((normal_model_one == JKR) || (normal_model_one == DMT))
|
||||
|
|
@ -894,10 +859,9 @@ 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;
|
||||
if (tangential_coeffs_one[0] == -1){
|
||||
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{
|
||||
} else {
|
||||
tangential_coeffs[i][j][0] = tangential_coeffs[j][i][0] = tangential_coeffs_one[0];
|
||||
}
|
||||
for (int k = 1; k < 3; k++)
|
||||
|
|
@ -949,23 +913,22 @@ void PairGranular::init_style()
|
|||
size_history = 3*tangential_history + 3*roll_history + twist_history;
|
||||
|
||||
//Determine location of tangential/roll/twist histories in array
|
||||
if (roll_history){
|
||||
if (roll_history) {
|
||||
if (tangential_history) roll_history_index = 3;
|
||||
else roll_history_index = 0;
|
||||
}
|
||||
if (twist_history){
|
||||
if (tangential_history){
|
||||
if (twist_history) {
|
||||
if (tangential_history) {
|
||||
if (roll_history) twist_history_index = 6;
|
||||
else twist_history_index = 3;
|
||||
}
|
||||
else{
|
||||
} else {
|
||||
if (roll_history) twist_history_index = 3;
|
||||
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){
|
||||
if (tangential_model[i][j] == TANGENTIAL_MINDLIN_RESCALE) {
|
||||
size_history += 1;
|
||||
roll_history_index += 1;
|
||||
twist_history_index += 1;
|
||||
|
|
@ -1049,24 +1012,23 @@ void PairGranular::init_style()
|
|||
int *type = atom->type;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
for (i = 0; i < nlocal; i++){
|
||||
for (i = 0; i < nlocal; i++) {
|
||||
double radius_cut = radius[i];
|
||||
if (mask[i] & freeze_group_bit){
|
||||
if (mask[i] & freeze_group_bit) {
|
||||
onerad_frozen[type[i]] = MAX(onerad_frozen[type[i]],radius_cut);
|
||||
}
|
||||
else{
|
||||
} else {
|
||||
onerad_dynamic[type[i]] = MAX(onerad_dynamic[type[i]],radius_cut);
|
||||
}
|
||||
}
|
||||
|
||||
MPI_Allreduce(&onerad_dynamic[1],&maxrad_dynamic[1],atom->ntypes,
|
||||
MPI_DOUBLE,MPI_MAX,world);
|
||||
MPI_DOUBLE,MPI_MAX,world);
|
||||
MPI_Allreduce(&onerad_frozen[1],&maxrad_frozen[1],atom->ntypes,
|
||||
MPI_DOUBLE,MPI_MAX,world);
|
||||
MPI_DOUBLE,MPI_MAX,world);
|
||||
|
||||
// set fix which stores history info
|
||||
|
||||
if (size_history > 0){
|
||||
if (size_history > 0) {
|
||||
int ifix = modify->find_fix("NEIGH_HISTORY");
|
||||
if (ifix < 0) error->all(FLERR,"Could not find pair fix neigh history ID");
|
||||
fix_history = (FixNeighHistory *) modify->fix[ifix];
|
||||
|
|
@ -1081,12 +1043,12 @@ double PairGranular::init_one(int i, int j)
|
|||
{
|
||||
double cutoff;
|
||||
|
||||
if (setflag[i][j] == 0){
|
||||
if (setflag[i][j] == 0) {
|
||||
if ((normal_model[i][i] != normal_model[j][j]) ||
|
||||
(damping_model[i][i] != damping_model[j][j]) ||
|
||||
(tangential_model[i][i] != tangential_model[j][j]) ||
|
||||
(roll_model[i][i] != roll_model[j][j]) ||
|
||||
(twist_model[i][i] != twist_model[j][j])){
|
||||
(twist_model[i][i] != twist_model[j][j])) {
|
||||
|
||||
char str[512];
|
||||
sprintf(str,"Granular pair style functional forms are different, cannot mix coefficients for types %d and %d. \nThis combination must be set explicitly via pair_coeff command.",i,j);
|
||||
|
|
@ -1105,12 +1067,12 @@ double PairGranular::init_one(int i, int j)
|
|||
for (int k = 0; k < 3; k++)
|
||||
tangential_coeffs[i][j][k] = tangential_coeffs[j][i][k] = mix_geom(tangential_coeffs[i][i][k], tangential_coeffs[j][j][k]);
|
||||
|
||||
if (roll_model[i][j] != ROLL_NONE){
|
||||
if (roll_model[i][j] != ROLL_NONE) {
|
||||
for (int k = 0; k < 3; k++)
|
||||
roll_coeffs[i][j][k] = roll_coeffs[j][i][k] = mix_geom(roll_coeffs[i][i][k], roll_coeffs[j][j][k]);
|
||||
}
|
||||
|
||||
if (twist_model[i][j] != TWIST_NONE && twist_model[i][j] != TWIST_MARSHALL){
|
||||
if (twist_model[i][j] != TWIST_NONE && twist_model[i][j] != TWIST_MARSHALL) {
|
||||
for (int k = 0; k < 3; k++)
|
||||
twist_coeffs[i][j][k] = twist_coeffs[j][i][k] = mix_geom(twist_coeffs[i][i][k], twist_coeffs[j][j][k]);
|
||||
}
|
||||
|
|
@ -1124,14 +1086,14 @@ double PairGranular::init_one(int i, int j)
|
|||
// we assign cutoff = max(cut[i][j]) for i,j such that cut[i][j] > 0.0.
|
||||
double pulloff;
|
||||
|
||||
if (cutoff_type[i][j] < 0 && cutoff_global < 0){
|
||||
if (cutoff_type[i][j] < 0 && cutoff_global < 0) {
|
||||
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){
|
||||
if (normal_model[i][j] == JKR) {
|
||||
pulloff = pulloff_distance(maxrad_dynamic[i], maxrad_dynamic[j], i, j);
|
||||
cutoff += pulloff;
|
||||
}
|
||||
|
|
@ -1143,8 +1105,7 @@ double PairGranular::init_one(int i, int j)
|
|||
if (normal_model[i][j] == JKR)
|
||||
pulloff = pulloff_distance(maxrad_dynamic[i], maxrad_frozen[j], i, j);
|
||||
cutoff = MAX(cutoff,maxrad_dynamic[i]+maxrad_frozen[j]+pulloff);
|
||||
}
|
||||
else { // radius info about either i or j does not exist (i.e. not present and not about to get poured; set to largest value to not interfere with neighbor list)
|
||||
} else { // radius info about either i or j does not exist (i.e. not present and not about to get poured; set to largest value to not interfere with neighbor list)
|
||||
double cutmax = 0.0;
|
||||
for (int k = 1; k <= atom->ntypes; k++) {
|
||||
cutmax = MAX(cutmax,2.0*maxrad_dynamic[k]);
|
||||
|
|
@ -1152,11 +1113,9 @@ double PairGranular::init_one(int i, int j)
|
|||
}
|
||||
cutoff = cutmax;
|
||||
}
|
||||
}
|
||||
else if (cutoff_type[i][j] > 0){
|
||||
} else if (cutoff_type[i][j] > 0) {
|
||||
cutoff = cutoff_type[i][j];
|
||||
}
|
||||
else if (cutoff_global > 0){
|
||||
} else if (cutoff_global > 0) {
|
||||
cutoff = cutoff_global;
|
||||
}
|
||||
|
||||
|
|
@ -1285,7 +1244,7 @@ double PairGranular::single(int i, int j, int itype, int jtype,
|
|||
|
||||
bool touchflag;
|
||||
E = normal_coeffs[itype][jtype][0];
|
||||
if (normal_model[itype][jtype] == JKR){
|
||||
if (normal_model[itype][jtype] == JKR) {
|
||||
E *= THREEQUARTERS;
|
||||
R2 = Reff*Reff;
|
||||
coh = normal_coeffs[itype][jtype][3];
|
||||
|
|
@ -1293,12 +1252,9 @@ double PairGranular::single(int i, int j, int itype, int jtype,
|
|||
delta_pulloff = a*a/Reff - 2*sqrt(M_PI*coh*a/E);
|
||||
dist_pulloff = radsum+delta_pulloff;
|
||||
touchflag = (rsq <= dist_pulloff*dist_pulloff);
|
||||
}
|
||||
else{
|
||||
touchflag = (rsq <= radsum*radsum);
|
||||
}
|
||||
} else 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;
|
||||
|
|
@ -1376,7 +1332,7 @@ double PairGranular::single(int i, int j, int itype, int jtype,
|
|||
|
||||
delta = radsum - r;
|
||||
dR = delta*Reff;
|
||||
if (normal_model[itype][jtype] == JKR){
|
||||
if (normal_model[itype][jtype] == JKR) {
|
||||
dR2 = dR*dR;
|
||||
t0 = coh*coh*R2*R2*E;
|
||||
t1 = PI27SQ*t0;
|
||||
|
|
@ -1392,12 +1348,11 @@ double PairGranular::single(int i, int j, int itype, int jtype,
|
|||
a2 = a*a;
|
||||
knfac = normal_coeffs[itype][jtype][0]*a;
|
||||
Fne = knfac*a2/Reff - TWOPI*a2*sqrt(4*coh*E/(M_PI*a));
|
||||
}
|
||||
else{
|
||||
} else {
|
||||
knfac = E;
|
||||
Fne = knfac*delta;
|
||||
a = sqrt(dR);
|
||||
if (normal_model[itype][jtype] != HOOKE){
|
||||
if (normal_model[itype][jtype] != HOOKE) {
|
||||
Fne *= a;
|
||||
knfac *= a;
|
||||
}
|
||||
|
|
@ -1405,13 +1360,11 @@ double PairGranular::single(int i, int j, int itype, int jtype,
|
|||
Fne -= 4*MY_PI*normal_coeffs[itype][jtype][3]*Reff;
|
||||
}
|
||||
|
||||
if (damping_model[itype][jtype] == VELOCITY){
|
||||
if (damping_model[itype][jtype] == VELOCITY) {
|
||||
damp_normal = normal_coeffs[itype][jtype][1];
|
||||
}
|
||||
else if (damping_model[itype][jtype] == VISCOELASTIC){
|
||||
} else if (damping_model[itype][jtype] == VISCOELASTIC) {
|
||||
damp_normal = normal_coeffs[itype][jtype][1]*a*meff;
|
||||
}
|
||||
else if (damping_model[itype][jtype] == TSUJI){
|
||||
} else if (damping_model[itype][jtype] == TSUJI) {
|
||||
damp_normal = normal_coeffs[itype][jtype][1]*sqrt(meff*knfac);
|
||||
}
|
||||
|
||||
|
|
@ -1423,7 +1376,7 @@ double PairGranular::single(int i, int j, int itype, int jtype,
|
|||
jnum = list->numneigh[i];
|
||||
jlist = list->firstneigh[i];
|
||||
|
||||
if (use_history){
|
||||
if (use_history) {
|
||||
allhistory = fix_history->firstvalue[i];
|
||||
for (int jj = 0; jj < jnum; jj++) {
|
||||
neighprev++;
|
||||
|
|
@ -1454,15 +1407,13 @@ double PairGranular::single(int i, int j, int itype, int jtype,
|
|||
vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3;
|
||||
vrel = sqrt(vrel);
|
||||
|
||||
if (normal_model[itype][jtype] == JKR){
|
||||
if (normal_model[itype][jtype] == JKR) {
|
||||
F_pulloff = 3*M_PI*coh*Reff;
|
||||
Fncrit = fabs(Fne + 2*F_pulloff);
|
||||
}
|
||||
else if (normal_model[itype][jtype] == DMT){
|
||||
} else if (normal_model[itype][jtype] == DMT) {
|
||||
F_pulloff = 4*M_PI*coh*Reff;
|
||||
Fncrit = fabs(Fne + 2*F_pulloff);
|
||||
}
|
||||
else{
|
||||
} else {
|
||||
Fncrit = fabs(Fntot);
|
||||
}
|
||||
|
||||
|
|
@ -1472,13 +1423,12 @@ double PairGranular::single(int i, int j, int itype, int jtype,
|
|||
k_tangential = tangential_coeffs[itype][jtype][0];
|
||||
damp_tangential = tangential_coeffs[itype][jtype][1]*damp_normal_prefactor;
|
||||
|
||||
if (tangential_history){
|
||||
if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN){
|
||||
if (tangential_history) {
|
||||
if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN) {
|
||||
k_tangential *= a;
|
||||
}
|
||||
else if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE){
|
||||
} else if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) {
|
||||
k_tangential *= a;
|
||||
if (a < history[3]){ //On unloading, rescale the shear displacements
|
||||
if (a < history[3]) { //On unloading, rescale the shear displacements
|
||||
double factor = a/history[3];
|
||||
history[0] *= factor;
|
||||
history[1] *= factor;
|
||||
|
|
@ -1507,8 +1457,7 @@ double PairGranular::single(int i, int j, int itype, int jtype,
|
|||
fs3 *= Fscrit/fs;
|
||||
} else fs1 = fs2 = fs3 = 0.0;
|
||||
}
|
||||
}
|
||||
else{ //Classic pair gran/hooke (no history)
|
||||
} 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;
|
||||
|
|
@ -1521,7 +1470,7 @@ double PairGranular::single(int i, int j, int itype, int jtype,
|
|||
// Rolling resistance
|
||||
//****************************************
|
||||
|
||||
if (roll_model[itype][jtype] != ROLL_NONE){
|
||||
if (roll_model[itype][jtype] != ROLL_NONE) {
|
||||
relrot1 = omega[i][0] - omega[j][0];
|
||||
relrot2 = omega[i][1] - omega[j][1];
|
||||
relrot3 = omega[i][2] - omega[j][2];
|
||||
|
|
@ -1565,14 +1514,13 @@ double PairGranular::single(int i, int j, int itype, int jtype,
|
|||
//****************************************
|
||||
// Twisting torque, including history effects
|
||||
//****************************************
|
||||
if (twist_model[itype][jtype] != TWIST_NONE){
|
||||
if (twist_model[itype][jtype] != TWIST_NONE) {
|
||||
magtwist = relrot1*nx + relrot2*ny + relrot3*nz; //Omega_T (eq 29 of Marshall)
|
||||
if (twist_model[itype][jtype] == TWIST_MARSHALL){
|
||||
if (twist_model[itype][jtype] == TWIST_MARSHALL) {
|
||||
k_twist = 0.5*k_tangential*a*a;; //eq 32
|
||||
damp_twist = 0.5*damp_tangential*a*a;
|
||||
mu_twist = TWOTHIRDS*a*tangential_coeffs[itype][jtype][2];;
|
||||
}
|
||||
else{
|
||||
} else {
|
||||
k_twist = twist_coeffs[itype][jtype][0];
|
||||
damp_twist = twist_coeffs[itype][jtype][1];
|
||||
mu_twist = twist_coeffs[itype][jtype][2];
|
||||
|
|
@ -1686,7 +1634,7 @@ double PairGranular::pulloff_distance(double radi, double radj, int itype, int j
|
|||
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){
|
||||
void PairGranular::transfer_history(double* source, double* target) {
|
||||
for (int i = 0; i < size_history; i++)
|
||||
target[i] = history_transfer_factors[i]*source[i];
|
||||
}
|
||||
|
|
|
|||
Loading…
Reference in New Issue