diff --git a/doc/src/fix_wall_gran.rst b/doc/src/fix_wall_gran.rst index 1eefc4aeae..57bec679e4 100644 --- a/doc/src/fix_wall_gran.rst +++ b/doc/src/fix_wall_gran.rst @@ -29,6 +29,7 @@ Syntax gamma_t = damping coefficient for collisions in tangential direction (1/time units or 1/time-distance units - see discussion below) xmu = static yield criterion (unitless value between 0.0 and 1.0e4) dampflag = 0 or 1 if tangential damping force is excluded or included + optional keyword = *limit_damping*, limit damping to prevent attractive interaction .. parsed-literal:: @@ -95,7 +96,8 @@ Specifically, delta = radius - r = overlap of particle with wall, m_eff = mass of particle, and the effective radius of contact = RiRj/Ri+Rj is set to the radius of the particle. -The parameters *Kn*\ , *Kt*\ , *gamma_n*, *gamma_t*, *xmu* and *dampflag* +The parameters *Kn*\ , *Kt*\ , *gamma_n*, *gamma_t*, *xmu*, *dampflag*, +and the optional keyword *limit_damping* have the same meaning and units as those specified with the :doc:`pair_style gran/\* ` commands. This means a NULL can be used for either *Kt* or *gamma_t* as described on that page. If a diff --git a/doc/src/fix_wall_gran_region.rst b/doc/src/fix_wall_gran_region.rst index 9ed5bb7bcc..0a252b162a 100644 --- a/doc/src/fix_wall_gran_region.rst +++ b/doc/src/fix_wall_gran_region.rst @@ -181,7 +181,8 @@ radius - r = overlap of particle with wall, m_eff = mass of particle, and the effective radius of contact is just the radius of the particle. -The parameters *Kn*\ , *Kt*\ , *gamma_n*, *gamma_t*, *xmu* and *dampflag* +The parameters *Kn*\ , *Kt*\ , *gamma_n*, *gamma_t*, *xmu*, *dampflag*, +and the optional keyword *limit_damping* have the same meaning and units as those specified with the :doc:`pair_style gran/\* ` commands. This means a NULL can be used for either *Kt* or *gamma_t* as described on that page. If a diff --git a/doc/src/pair_gran.rst b/doc/src/pair_gran.rst index d3f87839e8..fbcacb5c76 100644 --- a/doc/src/pair_gran.rst +++ b/doc/src/pair_gran.rst @@ -26,7 +26,7 @@ Syntax .. code-block:: LAMMPS - pair_style style Kn Kt gamma_n gamma_t xmu dampflag + pair_style style Kn Kt gamma_n gamma_t xmu dampflag keyword * style = *gran/hooke* or *gran/hooke/history* or *gran/hertz/history* * Kn = elastic constant for normal particle repulsion (force/distance units or pressure units - see discussion below) @@ -36,6 +36,13 @@ Syntax * xmu = static yield criterion (unitless value between 0.0 and 1.0e4) * dampflag = 0 or 1 if tangential damping force is excluded or included +* keyword = *limit_damping* + + .. parsed-literal:: + + *limit_damping* value = none + limit damping to prevent attractive interaction + .. note:: Versions of LAMMPS before 9Jan09 had different style names for @@ -54,6 +61,8 @@ Examples pair_style gran/hooke/history 200000.0 NULL 50.0 NULL 0.5 1 pair_style gran/hooke 200000.0 70000.0 50.0 30.0 0.5 0 + pair_style gran/hooke 200000.0 70000.0 50.0 30.0 0.5 0 limit_damping + Description """"""""""" @@ -208,6 +217,12 @@ potential is used as a sub-style of :doc:`pair_style hybrid `, then pair_coeff command to determine which atoms interact via a granular potential. +If two particles are moving away from each other while in contact, there +is a possibility that the particles could experience an effective attractive +force due to damping. If the *limit_damping* keyword is used, this option +will zero out the normal component of the force if there is an effective +attractive force. + ---------- .. include:: accel_styles.rst diff --git a/doc/src/pair_granular.rst b/doc/src/pair_granular.rst index 80663b7d49..432fd29584 100644 --- a/doc/src/pair_granular.rst +++ b/doc/src/pair_granular.rst @@ -24,7 +24,7 @@ Examples pair_coeff * * hooke 1000.0 50.0 tangential linear_history 500.0 1.0 0.4 damping mass_velocity pair_style granular - pair_coeff * * hertz 1000.0 50.0 tangential mindlin 1000.0 1.0 0.4 + pair_coeff * * hertz 1000.0 50.0 tangential mindlin 1000.0 1.0 0.4 limit_damping pair_style granular pair_coeff * * hertz/material 1e8 0.3 0.3 tangential mindlin_rescale NULL 1.0 0.4 damping tsuji @@ -623,6 +623,14 @@ Finally, the twisting torque on each particle is given by: ---------- +If two particles are moving away from each other while in contact, there +is a possibility that the particles could experience an effective attractive +force due to damping. If the optional *limit_damping* keyword is used, this option +will zero out the normal component of the force if there is an effective +attractive force. This keyword cannot be used with the JKR or DMT models. + +---------- + The *granular* pair style can reproduce the behavior of the *pair gran/\** styles with the appropriate settings (some very minor differences can be expected due to corrections in @@ -657,6 +665,12 @@ then LAMMPS will use that cutoff for the specified atom type combination, and automatically set pairwise cutoffs for the remaining atom types. +If two particles are moving away from each other while in contact, there +is a possibility that the particles could experience an effective attractive +force due to damping. If the *limit_damping* keyword is used, this option +will zero out the normal component of the force if there is an effective +attractive force. This keyword cannot be used with the JKR or DMT models. + ---------- .. include:: accel_styles.rst diff --git a/src/GRANULAR/fix_wall_gran.cpp b/src/GRANULAR/fix_wall_gran.cpp index 2efc189e6d..51dca15e7b 100644 --- a/src/GRANULAR/fix_wall_gran.cpp +++ b/src/GRANULAR/fix_wall_gran.cpp @@ -54,7 +54,8 @@ enum {XPLANE=0,YPLANE=1,ZPLANE=2,ZCYLINDER,REGION}; enum {NONE,CONSTANT,EQUAL}; enum {DAMPING_NONE, VELOCITY, MASS_VELOCITY, VISCOELASTIC, TSUJI}; enum {TANGENTIAL_NONE, TANGENTIAL_NOHISTORY, TANGENTIAL_HISTORY, - TANGENTIAL_MINDLIN, TANGENTIAL_MINDLIN_RESCALE}; + TANGENTIAL_MINDLIN, TANGENTIAL_MINDLIN_RESCALE, + TANGENTIAL_MINDLIN_FORCE, TANGENTIAL_MINDLIN_RESCALE_FORCE}; enum {TWIST_NONE, TWIST_SDS, TWIST_MARSHALL}; enum {ROLL_NONE, ROLL_SDS}; @@ -70,6 +71,7 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR,"Fix wall/gran requires atom style sphere"); create_attribute = 1; + limit_damping = 0; // set interaction style // disable bonded/history option for now @@ -90,7 +92,6 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : // wall/particle coefficients int iarg; - if (pairstyle != GRANULAR) { size_history = 3; if (narg < 11) error->all(FLERR,"Illegal fix wall/gran command"); @@ -118,6 +119,12 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : kt /= force->nktv2p; } iarg = 10; + + if (strcmp(arg[iarg],"limit_damping") == 0) { + limit_damping = 1; + iarg += 1; + } + } else { iarg = 4; damping_model = VISCOELASTIC; @@ -211,7 +218,9 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : iarg += 4; } 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) || + (strcmp(arg[iarg+1], "mindlin/force") == 0) || + (strcmp(arg[iarg+1], "mindlin_rescale/force") == 0)) { if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, " "not enough parameters provided for tangential model"); @@ -221,8 +230,14 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : tangential_model = TANGENTIAL_MINDLIN; else if (strcmp(arg[iarg+1], "mindlin_rescale") == 0) tangential_model = TANGENTIAL_MINDLIN_RESCALE; + else if (strcmp(arg[iarg+1], "mindlin/force") == 0) + tangential_model = TANGENTIAL_MINDLIN_FORCE; + else if (strcmp(arg[iarg+1], "mindlin_rescale/force") == 0) + tangential_model = TANGENTIAL_MINDLIN_RESCALE_FORCE; if ((tangential_model == TANGENTIAL_MINDLIN || - tangential_model == TANGENTIAL_MINDLIN_RESCALE) && + tangential_model == TANGENTIAL_MINDLIN_RESCALE || + tangential_model == TANGENTIAL_MINDLIN_FORCE || + tangential_model == TANGENTIAL_MINDLIN_RESCALE_FORCE) && (strcmp(arg[iarg+2], "NULL") == 0)) { if (normal_model == NORMAL_HERTZ || normal_model == NORMAL_HOOKE) { error->all(FLERR, "NULL setting for Mindlin tangential " @@ -293,13 +308,27 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : strcmp(arg[iarg], "zcylinder") == 0 || strcmp(arg[iarg], "region") == 0) { break; + } else if (strcmp(arg[iarg],"limit_damping") == 0) { + limit_damping = 1; + iarg += 1; } else { error->all(FLERR, "Illegal fix wall/gran command"); } } size_history = 3*tangential_history + 3*roll_history + twist_history; + //Unlike the pair style, the wall style does not have a 'touch' + //array. Hence, an additional entry in the history is used to + //determine if particles previously contacted for JKR cohesion purposes. if (normal_model == JKR) size_history += 1; - if (tangential_model == TANGENTIAL_MINDLIN_RESCALE) size_history += 1; + if (tangential_model == TANGENTIAL_MINDLIN_RESCALE || + tangential_model == TANGENTIAL_MINDLIN_RESCALE_FORCE) size_history += 1; + + if (limit_damping && normal_model == JKR) + error->all(FLERR,"Illegal pair_coeff command, " + "cannot limit damping with JRK model"); + if (limit_damping && normal_model == DMT) + error->all(FLERR,"Illegal pair_coeff command, " + "Cannot limit damping with DMT model"); } // wallstyle args @@ -475,36 +504,39 @@ void FixWallGran::init() 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; + if(pairstyle == GRANULAR) { + tangential_history_index = 0; + if (roll_history) { + if (tangential_history) roll_history_index = 3; + else roll_history_index = 0; } - else{ - if (roll_history) twist_history_index = 3; - else twist_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 (tangential_model == TANGENTIAL_MINDLIN_RESCALE || + tangential_model == TANGENTIAL_MINDLIN_RESCALE_FORCE) { + roll_history_index += 1; + twist_history_index += 1; } - } - if (normal_model == JKR) { - tangential_history_index += 1; - roll_history_index += 1; - twist_history_index += 1; - } - if (tangential_model == TANGENTIAL_MINDLIN_RESCALE) { - roll_history_index += 1; - twist_history_index += 1; - } - if (damping_model == TSUJI) { - double cor = normal_coeffs[1]; - 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); + 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); + } } } @@ -762,6 +794,7 @@ void FixWallGran::hooke(double rsq, double dx, double dy, double dz, damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radius-r)*rinv - damp; + if (limit_damping && (ccel < 0.0)) ccel = 0.0; // relative velocities @@ -854,6 +887,7 @@ void FixWallGran::hooke_history(double rsq, double dx, double dy, double dz, damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radius-r)*rinv - damp; + if (limit_damping && (ccel < 0.0)) ccel = 0.0; // relative velocities @@ -985,6 +1019,7 @@ void FixWallGran::hertz_history(double rsq, double dx, double dy, double dz, if (rwall == 0.0) polyhertz = sqrt((radius-r)*radius); else polyhertz = sqrt((radius-r)*radius*rwall/(rwall+radius)); ccel *= polyhertz; + if (limit_damping && (ccel < 0.0)) ccel = 0.0; // relative velocities @@ -1103,7 +1138,8 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz, double signtwist, magtwist, magtortwist, Mtcrit; double tortwist1, tortwist2, tortwist3; - double shrmag,rsht; + double shrmag,rsht,prjmag; + bool frameupdate; r = sqrt(rsq); E = normal_coeffs[0]; @@ -1178,11 +1214,18 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz, Fdamp = -damp_normal_prefactor*vnnr; Fntot = Fne + Fdamp; + if (limit_damping && (Fntot < 0.0)) Fntot = 0.0; //**************************************** // tangential force, including history effects //**************************************** + // For linear, mindlin, mindlin_rescale: + // history = cumulative tangential displacement + // + // For mindlin/force, mindlin_rescale/force: + // history = cumulative tangential elastic force + // tangential component vt1 = vr1 - vn1; vt2 = vr2 - vn2; @@ -1224,69 +1267,115 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz, int thist2 = thist1 + 1; if (tangential_history) { - if (tangential_model == TANGENTIAL_MINDLIN) { + if (tangential_model == TANGENTIAL_MINDLIN || + tangential_model == TANGENTIAL_MINDLIN_FORCE) { k_tangential *= a; } - else if (tangential_model == TANGENTIAL_MINDLIN_RESCALE) { + else if (tangential_model == + TANGENTIAL_MINDLIN_RESCALE || + tangential_model == + TANGENTIAL_MINDLIN_RESCALE_FORCE){ k_tangential *= a; - if (a < history[3]) { //On unloading, rescale the shear displacements + // on unloading, rescale the shear displacements/force + if (a < history[thist2+1]) { double factor = a/history[thist2+1]; history[thist0] *= factor; history[thist1] *= factor; history[thist2] *= factor; } } - shrmag = sqrt(history[thist0]*history[thist0] + - history[thist1]*history[thist1] + - history[thist2]*history[thist2]); + // 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) { - // if rhst == shrmag, contacting pair has rotated 90 deg in one step, - // in which case you deserve a crash! - scalefac = shrmag/(shrmag - rsht); + if (tangential_model == TANGENTIAL_MINDLIN_FORCE || + tangential_model == TANGENTIAL_MINDLIN_RESCALE_FORCE) + frameupdate = fabs(rsht) > EPSILON*Fscrit; + else + frameupdate = fabs(rsht)*k_tangential > EPSILON*Fscrit; + if (frameupdate) { + shrmag = sqrt(history[thist0]*history[thist0] + + history[thist1]*history[thist1] + + history[thist2]*history[thist2]); + // projection history[thist0] -= rsht*nx; history[thist1] -= rsht*ny; history[thist2] -= rsht*nz; + // also rescale to preserve magnitude + prjmag = sqrt(history[thist0]*history[thist0] + + history[thist1]*history[thist1] + history[thist2]*history[thist2]); + if (prjmag > 0) scalefac = shrmag/prjmag; + else scalefac = 0; history[thist0] *= scalefac; history[thist1] *= scalefac; history[thist2] *= scalefac; } // update history + if (tangential_model == TANGENTIAL_HISTORY || + tangential_model == TANGENTIAL_MINDLIN || + tangential_model == TANGENTIAL_MINDLIN_RESCALE) { history[thist0] += vtr1*dt; history[thist1] += vtr2*dt; history[thist2] += vtr3*dt; + } else{ + // tangential force + // see e.g. eq. 18 of Thornton et al, Pow. Tech. 2013, v223,p30-46 + history[thist0] -= k_tangential*vtr1*dt; + history[thist1] -= k_tangential*vtr2*dt; + history[thist2] -= k_tangential*vtr3*dt; + } + if (tangential_model == TANGENTIAL_MINDLIN_RESCALE || + tangential_model == TANGENTIAL_MINDLIN_RESCALE_FORCE) + history[thist2+1] = a; } // 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 (tangential_model == TANGENTIAL_HISTORY || + tangential_model == TANGENTIAL_MINDLIN || + tangential_model == TANGENTIAL_MINDLIN_RESCALE) { + fs1 = -k_tangential*history[thist0] - damp_tangential*vtr1; + fs2 = -k_tangential*history[thist1] - damp_tangential*vtr2; + fs3 = -k_tangential*history[thist2] - damp_tangential*vtr3; + } else { + fs1 = history[thist0] - damp_tangential*vtr1; + fs2 = history[thist1] - damp_tangential*vtr2; + fs3 = history[thist2] - damp_tangential*vtr3; + } // rescale frictional displacements and forces if needed Fscrit = tangential_coeffs[2] * Fncrit; fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); if (fs > Fscrit) { + shrmag = sqrt(history[thist0]*history[thist0] + + history[thist1]*history[thist1] + + history[thist2]*history[thist2]); 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); + if (tangential_model == TANGENTIAL_HISTORY || + tangential_model == TANGENTIAL_MINDLIN || + tangential_model == + TANGENTIAL_MINDLIN_RESCALE) { + 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); + } else { + history[thist0] = Fscrit*fs1/fs + damp_tangential*vtr1; + history[thist1] = Fscrit*fs2/fs + damp_tangential*vtr2; + history[thist2] = 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; + fs = damp_tangential*vrel; + if (vrel != 0.0) Ft = MIN(Fscrit,fs) / vrel; else Ft = 0.0; fs1 = -Ft*vtr1; fs2 = -Ft*vtr2; @@ -1297,14 +1386,14 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz, // rolling resistance //**************************************** - if (roll_model != ROLL_NONE || twist_model != NONE) { + if (roll_model != ROLL_NONE || twist_model != TWIST_NONE) { relrot1 = omega[0]; relrot2 = omega[1]; relrot3 = omega[2]; } if (roll_model != ROLL_NONE) { - - // rolling velocity, see eq. 31 of Wang et al, Particuology v 23, p 49 (2015) + // 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) @@ -1316,21 +1405,29 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz, int rhist1 = rhist0 + 1; int rhist2 = rhist1 + 1; - // rolling displacement - rollmag = sqrt(history[rhist0]*history[rhist0] + - history[rhist1]*history[rhist1] + - history[rhist2]*history[rhist2]); - - rolldotn = history[rhist0]*nx + history[rhist1]*ny + history[rhist2]*nz; + k_roll = roll_coeffs[0]; + damp_roll = roll_coeffs[1]; + Frcrit = roll_coeffs[2] * Fncrit; if (history_update) { - if (fabs(rolldotn) < EPSILON) rolldotn = 0; - if (rolldotn > 0) { // rotate into tangential plane - scalefac = rollmag/(rollmag - rolldotn); + rolldotn = history[rhist0]*nx + history[rhist1]*ny + history[rhist2]*nz; + frameupdate = fabs(rolldotn)*k_roll > EPSILON*Frcrit; + if (frameupdate) { // rotate into tangential plane + rollmag = sqrt(history[rhist0]*history[rhist0] + + history[rhist1]*history[rhist1] + + history[rhist2]*history[rhist2]); + // projection history[rhist0] -= rolldotn*nx; history[rhist1] -= rolldotn*ny; history[rhist2] -= rolldotn*nz; + // also rescale to preserve magnitude + prjmag = sqrt(history[rhist0]*history[rhist0] + + history[rhist1]*history[rhist1] + + history[rhist2]*history[rhist2]); + + if (prjmag > 0) scalefac = rollmag/prjmag; + else scalefac = 0; history[rhist0] *= scalefac; history[rhist1] *= scalefac; history[rhist2] *= scalefac; @@ -1340,17 +1437,16 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz, history[rhist2] += vrl3*dt; } - 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; // rescale frictional displacements and forces if needed - Frcrit = roll_coeffs[2] * Fncrit; - fr = sqrt(fr1*fr1 + fr2*fr2 + fr3*fr3); if (fr > Frcrit) { + rollmag = sqrt(history[rhist0]*history[rhist0] + + history[rhist1]*history[rhist1] + + history[rhist2]*history[rhist2]); if (rollmag != 0.0) { history[rhist0] = -1.0/k_roll*(Frcrit*fr1/fr + damp_roll*vrl1); history[rhist1] = -1.0/k_roll*(Frcrit*fr2/fr + damp_roll*vrl2); diff --git a/src/GRANULAR/fix_wall_gran.h b/src/GRANULAR/fix_wall_gran.h index bf4797e960..da48150fd0 100644 --- a/src/GRANULAR/fix_wall_gran.h +++ b/src/GRANULAR/fix_wall_gran.h @@ -73,6 +73,7 @@ class FixWallGran : public Fix { // for granular model choices int normal_model, damping_model; int tangential_model, roll_model, twist_model; + int limit_damping; // history flags int normal_history, tangential_history, roll_history, twist_history; diff --git a/src/GRANULAR/pair_gran_hertz_history.cpp b/src/GRANULAR/pair_gran_hertz_history.cpp index e7193140b9..56e9691fad 100644 --- a/src/GRANULAR/pair_gran_hertz_history.cpp +++ b/src/GRANULAR/pair_gran_hertz_history.cpp @@ -181,6 +181,7 @@ void PairGranHertzHistory::compute(int eflag, int vflag) ccel = kn*(radsum-r)*rinv - damp; polyhertz = sqrt((radsum-r)*radi*radj / radsum); ccel *= polyhertz; + if (limit_damping && (ccel < 0.0)) ccel = 0.0; // relative velocities @@ -277,7 +278,7 @@ void PairGranHertzHistory::compute(int eflag, int vflag) void PairGranHertzHistory::settings(int narg, char **arg) { - if (narg != 6) error->all(FLERR,"Illegal pair_style command"); + if (narg != 6 && narg != 7) error->all(FLERR,"Illegal pair_style command"); kn = utils::numeric(FLERR,arg[0],false,lmp); if (strcmp(arg[1],"NULL") == 0) kt = kn * 2.0/7.0; @@ -295,6 +296,12 @@ void PairGranHertzHistory::settings(int narg, char **arg) xmu < 0.0 || xmu > 10000.0 || dampflag < 0 || dampflag > 1) error->all(FLERR,"Illegal pair_style command"); + limit_damping = 0; + if (narg == 7) { + if (strcmp(arg[6], "limit_damping") == 0) limit_damping = 1; + else error->all(FLERR,"Illegal pair_style command"); + } + // convert Kn and Kt from pressure units to force/distance^2 kn /= force->nktv2p; @@ -388,6 +395,7 @@ double PairGranHertzHistory::single(int i, int j, int /*itype*/, int /*jtype*/, ccel = kn*(radsum-r)*rinv - damp; polyhertz = sqrt((radsum-r)*radi*radj / radsum); ccel *= polyhertz; + if (limit_damping && (ccel < 0.0)) ccel = 0.0; // relative velocities diff --git a/src/GRANULAR/pair_gran_hooke.cpp b/src/GRANULAR/pair_gran_hooke.cpp index 3875d12a0f..eb5ccec680 100644 --- a/src/GRANULAR/pair_gran_hooke.cpp +++ b/src/GRANULAR/pair_gran_hooke.cpp @@ -158,6 +158,7 @@ void PairGranHooke::compute(int eflag, int vflag) damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; + if (limit_damping && (ccel < 0.0)) ccel = 0.0; // relative velocities @@ -299,6 +300,7 @@ double PairGranHooke::single(int i, int j, int /*itype*/, int /*jtype*/, double damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; + if (limit_damping && (ccel < 0.0)) ccel = 0.0; // relative velocities diff --git a/src/GRANULAR/pair_gran_hooke_history.cpp b/src/GRANULAR/pair_gran_hooke_history.cpp index 4deb828d76..ac26ffbf1b 100644 --- a/src/GRANULAR/pair_gran_hooke_history.cpp +++ b/src/GRANULAR/pair_gran_hooke_history.cpp @@ -237,6 +237,7 @@ void PairGranHookeHistory::compute(int eflag, int vflag) damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; + if (limit_damping && (ccel < 0.0)) ccel = 0.0; // relative velocities @@ -356,7 +357,7 @@ void PairGranHookeHistory::allocate() void PairGranHookeHistory::settings(int narg, char **arg) { - if (narg != 6) error->all(FLERR,"Illegal pair_style command"); + if (narg != 6 && narg != 7) error->all(FLERR,"Illegal pair_style command"); kn = utils::numeric(FLERR,arg[0],false,lmp); if (strcmp(arg[1],"NULL") == 0) kt = kn * 2.0/7.0; @@ -370,6 +371,12 @@ void PairGranHookeHistory::settings(int narg, char **arg) dampflag = utils::inumeric(FLERR,arg[5],false,lmp); if (dampflag == 0) gammat = 0.0; + limit_damping = 0; + if (narg == 7) { + if (strcmp(arg[6], "limit_damping") == 0) limit_damping = 1; + else error->all(FLERR,"Illegal pair_style command"); + } + if (kn < 0.0 || kt < 0.0 || gamman < 0.0 || gammat < 0.0 || xmu < 0.0 || xmu > 10000.0 || dampflag < 0 || dampflag > 1) error->all(FLERR,"Illegal pair_style command"); @@ -686,6 +693,7 @@ double PairGranHookeHistory::single(int i, int j, int /*itype*/, int /*jtype*/, damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; + if(limit_damping && (ccel < 0.0)) ccel = 0.0; // relative velocities diff --git a/src/GRANULAR/pair_gran_hooke_history.h b/src/GRANULAR/pair_gran_hooke_history.h index c27ce8a9af..f511d30237 100644 --- a/src/GRANULAR/pair_gran_hooke_history.h +++ b/src/GRANULAR/pair_gran_hooke_history.h @@ -49,6 +49,7 @@ class PairGranHookeHistory : public Pair { double dt; int freeze_group_bit; int history; + int limit_damping; int neighprev; double *onerad_dynamic,*onerad_frozen; diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index bcf262aa2c..52971f920f 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -81,6 +81,8 @@ PairGranular::PairGranular(LAMMPS *lmp) : Pair(lmp) maxrad_dynamic = nullptr; maxrad_frozen = nullptr; + limit_damping = nullptr; + history_transfer_factors = nullptr; dt = update->dt; @@ -130,6 +132,7 @@ PairGranular::~PairGranular() memory->destroy(tangential_model); memory->destroy(roll_model); memory->destroy(twist_model); + memory->destroy(limit_damping); delete [] onerad_dynamic; delete [] onerad_frozen; @@ -364,12 +367,13 @@ void PairGranular::compute(int eflag, int vflag) damp_normal = a*meff; } else if (damping_model[itype][jtype] == TSUJI) { damp_normal = sqrt(meff*knfac); - } + } else damp_normal = 0.0; damp_normal_prefactor = normal_coeffs[itype][jtype][1]*damp_normal; Fdamp = -damp_normal_prefactor*vnnr; Fntot = Fne + Fdamp; + if (limit_damping[itype][jtype] && (Fntot < 0.0)) Fntot = 0.0; //**************************************** // tangential force, including history effects @@ -540,20 +544,17 @@ void PairGranular::compute(int eflag, int vflag) relrot1 = omega[i][0] - omega[j][0]; relrot2 = omega[i][1] - omega[j][1]; relrot3 = omega[i][2] - omega[j][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) - // - 0.5*((radj-radi)/radsum)*vtr1; - // - 0.5*((radj-radi)/radsum)*vtr2; - // - 0.5*((radj-radi)/radsum)*vtr3; } //**************************************** // rolling resistance //**************************************** if (roll_model[itype][jtype] != ROLL_NONE) { + // 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); vrl2 = Reff*(relrot3*nx - relrot1*nz); vrl3 = Reff*(relrot1*ny - relrot2*nx); @@ -740,6 +741,7 @@ void PairGranular::allocate() memory->create(tangential_model,n+1,n+1,"pair:tangential_model"); memory->create(roll_model,n+1,n+1,"pair:roll_model"); memory->create(twist_model,n+1,n+1,"pair:twist_model"); + memory->create(limit_damping,n+1,n+1,"pair:limit_damping"); onerad_dynamic = new double[n+1]; onerad_frozen = new double[n+1]; @@ -793,6 +795,7 @@ void PairGranular::coeff(int narg, char **arg) roll_model_one = ROLL_NONE; twist_model_one = TWIST_NONE; damping_model_one = VISCOELASTIC; + int ld_flag = 0; int iarg = 2; while (iarg < narg) { @@ -967,12 +970,15 @@ void PairGranular::coeff(int narg, char **arg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters"); cutoff_one = utils::numeric(FLERR,arg[iarg+1],false,lmp); iarg += 2; - } else error->all(FLERR, "Illegal pair coeff command"); + } else if (strcmp(arg[iarg], "limit_damping") == 0) { + ld_flag = 1; + iarg += 1; + } else error->all(FLERR, "Illegal pair_coeff command"); } // error not to specify normal or tangential model if ((normal_model_one < 0) || (tangential_model_one < 0)) - error->all(FLERR, "Illegal pair coeff command, " + error->all(FLERR, "Illegal pair_coeff command, " "must specify normal or tangential contact model"); int count = 0; @@ -984,6 +990,14 @@ void PairGranular::coeff(int narg, char **arg) 27.467*powint(cor,4)-18.022*powint(cor,5)+4.8218*powint(cor,6); } else damp = normal_coeffs_one[1]; + if (ld_flag && normal_model_one == JKR) + error->all(FLERR,"Illegal pair_coeff command, " + "Cannot limit damping with JKR model"); + + if (ld_flag && normal_model_one == DMT) + error->all(FLERR,"Illegal pair_coeff command, " + "Cannot limit damping with DMT model"); + 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; @@ -1026,11 +1040,14 @@ void PairGranular::coeff(int narg, char **arg) cutoff_type[i][j] = cutoff_type[j][i] = cutoff_one; + limit_damping[i][j] = limit_damping[j][i] = ld_flag; + setflag[i][j] = 1; count++; } } + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); } @@ -1303,6 +1320,7 @@ void PairGranular::write_restart(FILE *fp) fwrite(&tangential_model[i][j],sizeof(int),1,fp); fwrite(&roll_model[i][j],sizeof(int),1,fp); fwrite(&twist_model[i][j],sizeof(int),1,fp); + fwrite(&limit_damping[i][j],sizeof(int),1,fp); fwrite(normal_coeffs[i][j],sizeof(double),4,fp); fwrite(tangential_coeffs[i][j],sizeof(double),3,fp); fwrite(roll_coeffs[i][j],sizeof(double),3,fp); @@ -1333,6 +1351,7 @@ void PairGranular::read_restart(FILE *fp) utils::sfread(FLERR,&tangential_model[i][j],sizeof(int),1,fp,nullptr,error); utils::sfread(FLERR,&roll_model[i][j],sizeof(int),1,fp,nullptr,error); utils::sfread(FLERR,&twist_model[i][j],sizeof(int),1,fp,nullptr,error); + utils::sfread(FLERR,&limit_damping[i][j],sizeof(int),1,fp,nullptr,error); utils::sfread(FLERR,normal_coeffs[i][j],sizeof(double),4,fp,nullptr,error); utils::sfread(FLERR,tangential_coeffs[i][j],sizeof(double),3,fp,nullptr,error); utils::sfread(FLERR,roll_coeffs[i][j],sizeof(double),3,fp,nullptr,error); @@ -1344,6 +1363,7 @@ void PairGranular::read_restart(FILE *fp) MPI_Bcast(&tangential_model[i][j],1,MPI_INT,0,world); MPI_Bcast(&roll_model[i][j],1,MPI_INT,0,world); MPI_Bcast(&twist_model[i][j],1,MPI_INT,0,world); + MPI_Bcast(&limit_damping[i][j],1,MPI_INT,0,world); MPI_Bcast(normal_coeffs[i][j],4,MPI_DOUBLE,0,world); MPI_Bcast(tangential_coeffs[i][j],3,MPI_DOUBLE,0,world); MPI_Bcast(roll_coeffs[i][j],3,MPI_DOUBLE,0,world); @@ -1529,6 +1549,7 @@ double PairGranular::single(int i, int j, int itype, int jtype, Fdamp = -damp_normal_prefactor*vnnr; Fntot = Fne + Fdamp; + if (limit_damping[itype][jtype] && (Fntot < 0.0)) Fntot = 0.0; jnum = list->numneigh[i]; jlist = list->firstneigh[i]; diff --git a/src/GRANULAR/pair_granular.h b/src/GRANULAR/pair_granular.h index 7ef4f2af98..49eb457f0d 100644 --- a/src/GRANULAR/pair_granular.h +++ b/src/GRANULAR/pair_granular.h @@ -70,6 +70,7 @@ class PairGranular : public Pair { // model choices int **normal_model, **damping_model; int **tangential_model, **roll_model, **twist_model; + int **limit_damping; // history flags int normal_history, tangential_history, roll_history, twist_history; diff --git a/src/KOKKOS/pair_gran_hooke_history_kokkos.cpp b/src/KOKKOS/pair_gran_hooke_history_kokkos.cpp index 81c82ab826..01ac1ea54f 100644 --- a/src/KOKKOS/pair_gran_hooke_history_kokkos.cpp +++ b/src/KOKKOS/pair_gran_hooke_history_kokkos.cpp @@ -392,6 +392,7 @@ void PairGranHookeHistoryKokkos::operator()(TagPairGranHookeHistoryC F_FLOAT damp = meff*gamman*vnnr*rsqinv; F_FLOAT ccel = kn*(radsum-r)*rinv - damp; + if(limit_damping && (ccel < 0.0)) ccel = 0.0; // relative velocities diff --git a/src/USER-OMP/pair_gran_hertz_history_omp.cpp b/src/USER-OMP/pair_gran_hertz_history_omp.cpp index 2187dcbe90..200075a316 100644 --- a/src/USER-OMP/pair_gran_hertz_history_omp.cpp +++ b/src/USER-OMP/pair_gran_hertz_history_omp.cpp @@ -224,6 +224,7 @@ void PairGranHertzHistoryOMP::eval(int iifrom, int iito, ThrData * const thr) ccel = kn*(radsum-r)*rinv - damp; polyhertz = sqrt((radsum-r)*radi*radj / radsum); ccel *= polyhertz; + if (limit_damping && (ccel < 0.0)) ccel = 0.0; // relative velocities diff --git a/src/USER-OMP/pair_gran_hooke_history_omp.cpp b/src/USER-OMP/pair_gran_hooke_history_omp.cpp index 7894d1e963..0ba9607bd9 100644 --- a/src/USER-OMP/pair_gran_hooke_history_omp.cpp +++ b/src/USER-OMP/pair_gran_hooke_history_omp.cpp @@ -223,6 +223,7 @@ void PairGranHookeHistoryOMP::eval(int iifrom, int iito, ThrData * const thr) damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; + if (limit_damping && (ccel < 0.0)) ccel = 0.0; // relative velocities diff --git a/src/USER-OMP/pair_gran_hooke_omp.cpp b/src/USER-OMP/pair_gran_hooke_omp.cpp index 08ae07c6f5..451b074fa9 100644 --- a/src/USER-OMP/pair_gran_hooke_omp.cpp +++ b/src/USER-OMP/pair_gran_hooke_omp.cpp @@ -197,6 +197,7 @@ void PairGranHookeOMP::eval(int iifrom, int iito, ThrData * const thr) damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; + if (limit_damping && (ccel < 0.0)) ccel = 0.0; // relative velocities