diff --git a/doc/src/dump.rst b/doc/src/dump.rst index 4d272c940b..4406c827bf 100644 --- a/doc/src/dump.rst +++ b/doc/src/dump.rst @@ -81,6 +81,7 @@ Syntax q, mux, muy, muz, mu, radius, diameter, omegax, omegay, omegaz, angmomx, angmomy, angmomz, tqx, tqy, tqz, + heatflux, temperature, c_ID, c_ID[I], f_ID, f_ID[I], v_name, i_name, d_name, i2_name[I], d2_name[I] @@ -103,10 +104,12 @@ Syntax q = atom charge mux,muy,muz = orientation of dipole moment of atom mu = magnitude of dipole moment of atom - radius,diameter = radius,diameter of spherical particle + radius,diameter = radius, diameter of spherical particle omegax,omegay,omegaz = angular velocity of spherical particle angmomx,angmomy,angmomz = angular momentum of aspherical particle tqx,tqy,tqz = torque on finite-size particles + heatflux = flux of heat into particle + temperature = temperature of particle c_ID = per-atom vector calculated by a compute with ID c_ID[I] = Ith column of per-atom array calculated by a compute with ID, I can include wildcard (see below) f_ID = per-atom vector calculated by a fix with ID diff --git a/doc/src/pair_granular.rst b/doc/src/pair_granular.rst index 6f84b0d9c7..9a99558653 100644 --- a/doc/src/pair_granular.rst +++ b/doc/src/pair_granular.rst @@ -37,6 +37,9 @@ Examples pair_coeff 1 1 dmt 1000.0 50.0 0.3 0.0 tangential mindlin NULL 0.5 0.5 rolling sds 500.0 200.0 0.5 twisting marshall pair_coeff 2 2 dmt 1000.0 50.0 0.3 10.0 tangential mindlin NULL 0.5 0.1 rolling sds 500.0 200.0 0.1 twisting marshall + pair_style granular + pair_coeff * * hertz 1000.0 50.0 tangential mindlin 1000.0 1.0 0.4 heat 0.1 + Description """"""""""" @@ -631,6 +634,13 @@ attractive force. This keyword cannot be used with the JKR or DMT models. ---------- +The optional *heat* keyword enables heat conduction and it must be followed by +a non-negative numeric value for the conductivity. Note that the *heat* keyword +must be included in either all or none of of the *pair_coeff* calls. See +:doc:`fix temp/integrate ` for more information on this option. + +---------- + 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 diff --git a/src/GRANULAR/atom_vec_sphere_temp.cpp b/src/GRANULAR/atom_vec_sphere_temp.cpp index d4e30e3154..a455a9ddc8 100644 --- a/src/GRANULAR/atom_vec_sphere_temp.cpp +++ b/src/GRANULAR/atom_vec_sphere_temp.cpp @@ -28,6 +28,7 @@ using namespace MathConst; AtomVecSphereTemp::AtomVecSphereTemp(LAMMPS *lmp) : AtomVec(lmp) { mass_type = PER_ATOM; + forceclearflag = 1; molecular = Atom::ATOMIC; atom->sphere_flag = 1; @@ -107,6 +108,16 @@ void AtomVecSphereTemp::grow_pointers() omega = atom->omega; } +/* ---------------------------------------------------------------------- + clear extra forces starting at atom N + nbytes = # of bytes to clear for a per-atom vector +------------------------------------------------------------------------- */ + +void AtomVecSphereTemp::force_clear(int n, size_t nbytes) +{ + memset(&heatflux[n], 0, nbytes); +} + /* ---------------------------------------------------------------------- initialize non-zero atom quantities ------------------------------------------------------------------------- */ diff --git a/src/GRANULAR/atom_vec_sphere_temp.h b/src/GRANULAR/atom_vec_sphere_temp.h index f42b844886..70e047a8cc 100644 --- a/src/GRANULAR/atom_vec_sphere_temp.h +++ b/src/GRANULAR/atom_vec_sphere_temp.h @@ -31,6 +31,7 @@ class AtomVecSphereTemp : public AtomVec { void init() override; void grow_pointers() override; + void force_clear(int, size_t) override; void create_atom_post(int) override; void data_atom_post(int) override; void pack_data_pre(int) override; diff --git a/src/GRANULAR/contact.cpp b/src/GRANULAR/contact.cpp index 3cea933370..b4b335a1aa 100644 --- a/src/GRANULAR/contact.cpp +++ b/src/GRANULAR/contact.cpp @@ -17,28 +17,18 @@ and torques based on contact geometry */ -#include #include "contact.h" +#include "math_const.h" +#include "math_extra.h" +#include "pointers.h" -namespace LAMMPS_NS { - namespace Contact{ +#include -enum {HOOKE, HERTZ, HERTZ_MATERIAL, DMT, JKR}; -enum {VELOCITY, MASS_VELOCITY, VISCOELASTIC, TSUJI}; -enum {TANGENTIAL_NOHISTORY, TANGENTIAL_HISTORY, - TANGENTIAL_MINDLIN, TANGENTIAL_MINDLIN_RESCALE, - TANGENTIAL_MINDLIN_FORCE, TANGENTIAL_MINDLIN_RESCALE_FORCE}; -enum {TWIST_NONE, TWIST_SDS, TWIST_MARSHALL}; -enum {ROLL_NONE, ROLL_SDS}; +using namespace LAMMPS_NS; +using namespace MathExtra; +using namespace MathConst; -#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 (4.0/3.0) // 4/3 -#define THREEQUARTERS 0.75 // 3/4 - -#define EPSILON 1e-10 +namespace Contact { ContactModel::ContactModel() { @@ -69,8 +59,8 @@ bool ContactModel::check_contact() { check_flag = 1; - MathExtra::sub3(xi, xj, dx); - rsq = MathExtra::lensq3(dx); + sub3(xi, xj, dx); + rsq = lensq3(dx); radsum = radi + radj; Reff = radi*radj/radsum; @@ -78,7 +68,7 @@ bool ContactModel::check_contact() if (normal_model == JKR) touch = touch_JKR(touch); else touch = (rsq < radsum*radsum); - return touch + return touch; } /* ---------------------------------------------------------------------- */ @@ -91,57 +81,57 @@ void ContactModel::prep_contact() if (check_flag != 1) touch = check_contact(); if (!touch) return; + double temp[3]; + // Set flags mindlin_rescale = mindlin_force = 0; if (tangential_model == TANGENTIAL_MINDLIN_RESCALE || tangential_model == TANGENTIAL_MINDLIN_RESCALE_FORCE) - mindlin_rescale = 1 + mindlin_rescale = 1; if (tangential_model == TANGENTIAL_MINDLIN_FORCE || tangential_model == TANGENTIAL_MINDLIN_RESCALE_FORCE) - mindlin_force = 1 - - double temp[3]; + mindlin_force = 1; // Standard geometric quantities r = sqrt(rsq); rinv = 1.0/r; delta = radsum - r; - dR = delta*Reff - MathExtra::scale3(rinv, dx, nx); + dR = delta*Reff; + scale3(rinv, dx, nx); // relative translational velocity - MathExtra::sub3(v[i], v[j], vr); + sub3(vi, vj, vr); // normal component - vnnr = MathExtra::dot3(vr, nx); //v_R . n - MathExtra::scale3(vnnr, nx, vn); + vnnr = dot3(vr, nx); //v_R . n + scale3(vnnr, nx, vn); // tangential component - MathExtra::sub3(vr, vn, vt); + sub3(vr, vn, vt); // relative rotational velocity - MathExtra::scaleadd3(radi, omegai, radj, omegaj, wr); + scaleadd3(radi, omegai, radj, omegaj, wr); // relative tangential velocities - MathExtra::cross3(wr, nx, temp); - MathExtra::sub3(vt, temp, vtr); - vrel = MathExtra::len(vtr); + cross3(wr, nx, temp); + sub3(vt, temp, vtr); + vrel = len3(vtr); - if (roll_model != NONE || twist_model != NONE) - MathExtra::sub3(omega[i], omega[j], relrot); + if (roll_model != ROLL_NONE || twist_model != TWIST_NONE) + sub3(omegai, omegaj, relrot); - if (roll_model != NONE) { + if (roll_model != 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) - vrl[0] = Reff * (relrot[1] * n[2] - relrot[2] * nx[1]); - vrl[1] = Reff * (relrot[2] * n[0] - relrot[0] * nx[2]); - vrl[2] = Reff * (relrot[0] * n[1] - relrot[1] * nx[0]); + vrl[0] = Reff * (relrot[1] * nx[2] - relrot[2] * nx[1]); + vrl[1] = Reff * (relrot[2] * nx[0] - relrot[0] * nx[2]); + vrl[2] = Reff * (relrot[0] * nx[1] - relrot[1] * nx[0]); } - if (twist_model != NONE) { + if (twist_model != TWIST_NONE) { // omega_T (eq 29 of Marshall) - magtwist = MathExtra::dot3(relrot, nx); + magtwist = dot3(relrot, nx); } } @@ -154,23 +144,23 @@ void ContactModel::calculate_forces(double *forces, double *torquesi, double *to if (prep_flag != 1) prep_contact(); if (!touch) { forces[0] = forces[1] = forces[2] = 0.0; - return + return; } //********************************************** // normal forces //********************************************** - // Also calculates: a, knfac, Fncrit + // Also calculates: a, knfac, Fncrit (for JKR or DMT) double Fne; if (normal_model == JKR) { Fne = normal_JKR(); } else if (normal_model == DMT) { Fne = normal_DMT(); } else if (normal_model == HERTZ || normal_model == HERTZ_MATERIAL) { - Fne = normal_hertz(); + Fne = normal_Hertz(); } else { - Fne = normal_hooke(); + Fne = normal_Hooke(); } // NOTE: consider restricting Hooke to only have @@ -178,9 +168,11 @@ void ContactModel::calculate_forces(double *forces, double *torquesi, double *to // Also calculates: damp_normal_prefactor double Fdamp = normal_damping(); - double Fntot = Fne + Fdamp; + Fntot = Fne + Fdamp; if (limit_damping && (Fntot < 0.0)) Fntot = 0.0; + if (normal_model != JKR && normal_model != DMT) Fncrit = fabs(Fntot); + //********************************************** // tangential force, including history effects //********************************************** @@ -190,15 +182,18 @@ void ContactModel::calculate_forces(double *forces, double *torquesi, double *to // For mindlin/force, mindlin_rescale/force: // history = cumulative tangential elastic force - Fncrit = critical_normal(Fne, Fntot, geom, model); Fscrit = mu_tang * Fncrit; - if (tangential_model == TANGENTIAL_NOHISTORY) { - tangential_no_history(); + if (tangential_model == TANGENTIAL_MINDLIN || + tangential_model == TANGENTIAL_MINDLIN_FORCE || + tangential_model == TANGENTIAL_MINDLIN_RESCALE || + tangential_model == TANGENTIAL_MINDLIN_RESCALE_FORCE) { + tangential_mindlin(history); } else if (tangential_model == TANGENTIAL_HISTORY) { tangential_history(history); + } else { + tangential_no_history(); } - tangential_forces(history); //********************************************** // rolling force @@ -212,50 +207,63 @@ void ContactModel::calculate_forces(double *forces, double *torquesi, double *to //**************************************** if (twist_model == TWIST_MARSHALL) { - twist_marshall(history); + twisting_marshall(history); } else if (twist_model == TWIST_SDS) { - twist_SDS(history); + twisting_SDS(history); } //********************************************** // sum contributions //********************************************** - MathExtra::scale3(Fntot, nx, forces); - MathExtra::add3(forces, fs, forces); + scale3(Fntot, nx, forces); + add3(forces, fs, forces); - MathExtra::cross3(nx, fs, torquesi); - MathExtra::copy3(torquesi, torquesj); + cross3(nx, fs, torquesi); + copy3(torquesi, torquesj); double dist_to_contact = radi-0.5*delta; - MathExtra::scale3(dist_to_contact, torquesi); + scale3(dist_to_contact, torquesi); dist_to_contact = radj-0.5*delta; - MathExtra::scale3(dist_to_contact, torquesj); + scale3(dist_to_contact, torquesj); double torroll[3]; if (roll_model != ROLL_NONE) { - MathExtra::cross3(nx, fr, torroll); - MathExtra::scale3(Reff, torroll); - MathExtra::add3(torquesi, torroll, torquesi); - MathExtra::sub3(torquesj, torroll, torquesj); + cross3(nx, fr, torroll); + scale3(Reff, torroll); + add3(torquesi, torroll, torquesi); + sub3(torquesj, torroll, torquesj); } double tortwist[3]; - if (twist_model != NONE_TWIST) { - MathExtra::scale3(magtortwist, nx, tortwist); - MathExtra::add3(torquesi, tortwist, torquesi) - MathExtra::sub3(torquesj, tortwist, torquesj) + if (twist_model != TWIST_NONE) { + scale3(magtortwist, nx, tortwist); + add3(torquesi, tortwist, torquesi); + sub3(torquesj, tortwist, torquesj); } } /* ---------------------------------------------------------------------- */ -double ContactModel::touch_JKR(int touch) +double ContactModel::calculate_heat() +{ + double dT = Ti - Tj; + double Hc; + + //Dan is Emod the youngs modulus for all models? or do I need to scale? + Hc = 2 * conductivity * pow(FOURTHIRDS * Fntot * Reff / Emod, ONETHIRD); + + return Hc * dT; +} + +/* ---------------------------------------------------------------------- */ + +bool ContactModel::touch_JKR(int touch) { double Escaled, R2, delta_pulloff, dist_pulloff; bool touchflag; - Escaled = E*THREEQUARTERS; + Escaled = k_norm * THREEQUARTERS; if (touch) { R2 = Reff * Reff; a = cbrt(9.0 * MY_PI * cohesion * R2 / (4 * Escaled)); @@ -272,66 +280,68 @@ double ContactModel::touch_JKR(int touch) double ContactModel::normal_JKR() { - double R2, dR2, t0, t1, t2, t3, t4, t5, t6; - double sqrt1, sqrt2, sqrt3, a, a2, F_pulloff; + double Escaled, R2, dR2, t0, t1, t2, t3, t4, t5, t6; + double sqrt1, sqrt2, sqrt3, a, a2, F_pulloff, Fne; + + Escaled = k_norm * THREEQUARTERS; R2 = Reff * Reff; dR2 = dR * dR; - t0 = cohesion * cohesion * R2 * R2 * E; + t0 = cohesion * cohesion * R2 * R2 * Escaled; t1 = PI27SQ*t0; - t2 = 8 * dR * dR2 * E * E * E; - t3 = 4 * dR2 * E; + t2 = 8 * dR * dR2 * Escaled * Escaled * Escaled; + t3 = 4 * dR2 * Escaled; // in case sqrt(0) < 0 due to precision issues sqrt1 = MAX(0, t0 * (t1 + 2 * t2)); t4 = cbrt(t1 + t2 + THREEROOT3 * MY_PI * sqrt(sqrt1)); - t5 = t3 / t4 + t4 / E; + t5 = t3 / t4 + t4 / Escaled; sqrt2 = MAX(0, 2 * dR + t5); t6 = sqrt(sqrt2); - sqrt3 = MAX(0, 4 * dR - t5 + SIXROOT6 * cohesion * MY_PI * R2 / (E * t6)); + sqrt3 = MAX(0, 4 * dR - t5 + SIXROOT6 * cohesion * MY_PI * R2 / (Escaled * t6)); a = INVROOT6 * (t6 + sqrt(sqrt3)); a2 = a * a; - double Fne = E * a * a2 / Reff - MY_2PI * a2 * sqrt(4 * cohesion * E / (MY_PI * a)); - double F_pulloff = 3 * MY_PI * cohesion * Reff; + Fne = Escaled * a * a2 / Reff - MY_2PI * a2 * sqrt(4 * cohesion * Escaled / (MY_PI * a)); + F_pulloff = 3 * MY_PI * cohesion * Reff; - knfac = E * a; - Fncrit = fabs(Fne + 2*F_pulloff); + knfac = Escaled * a; + Fncrit = fabs(Fne + 2 * F_pulloff); return Fne; } /* ---------------------------------------------------------------------- */ -double ContactModel::normal_DMT(double &Fne) +double ContactModel::normal_DMT() { a = sqrt(dR); - double Fne = a * E * delta; + double Fne = a * k_norm * delta; Fne -= 4 * MY_PI * cohesion * Reff; - F_pulloff = 4 * MY_PI * cohesion * Reff; + double F_pulloff = 4 * MY_PI * cohesion * Reff; - knfac = E * a; - Fncrit = fabs(Fne + 2*F_pulloff); + knfac = k_norm * a; + Fncrit = fabs(Fne + 2 * F_pulloff); return Fne; } /* ---------------------------------------------------------------------- */ -double ContactModel::normal_Hertz(double &Fne) +double ContactModel::normal_Hertz() { a = sqrt(dR); - double Fne = E * delta * a; + double Fne = k_norm * delta * a; - knfac = E * a; + knfac = k_norm * a; return Fne; } /* ---------------------------------------------------------------------- */ -double ContactModel::normal_Hooke(double &Fne) +double ContactModel::normal_Hooke() { a = sqrt(dR); - double Fne = E * delta; + double Fne = k_norm * delta; - knfac = E; + knfac = k_norm; return Fne; } @@ -356,7 +366,7 @@ double ContactModel::normal_damping() /* ---------------------------------------------------------------------- */ -double ContactModel::tangential_no_history() +void ContactModel::tangential_no_history() { double gamma_scaled = gamma_tang * damp_normal_prefactor; double fsmag, Ft; @@ -367,100 +377,100 @@ double ContactModel::tangential_no_history() else Ft = 0.0; Ft = -Ft; - MathExtra::scale3(Ft, vtr, fs); + scale3(Ft, vtr, fs); } /* ---------------------------------------------------------------------- */ -double ContactModel::tangential_forces(double *history) +void ContactModel::tangential_history(double *history) { double gamma_scaled = gamma_tang * damp_normal_prefactor; double k = k_tang; int frame_update = 0; - double fsmag, rsht, shrmag, prjmag, temp_dbl, temp_array[3]; + double magfs, rsht, shrmag, prjmag, temp_dbl, temp_array[3]; // rotate and update displacements / force. // see e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235 if (history_update) { - rsht = MathExtra::dot(history, nx); - frameupdate = fabs(rsht) * k > EPSILON * Fscrit; + rsht = dot3(history, nx); + frame_update = fabs(rsht) * k > EPSILON * Fscrit; - if (frameupdate) { - shrmag = MathExtra::len3(history); + if (frame_update) { + shrmag = len3(history); // projection - MathExtra::scale3(rsht, nx, history); + scale3(rsht, nx, history); // also rescale to preserve magnitude - prjmag = MathExtra::len3(history); - if (prjmag > 0) temp_double = shrmag / prjmag; - else temp_double = 0; - MathExtra::scale3(temp_double, history); + prjmag = len3(history); + if (prjmag > 0) temp_dbl = shrmag / prjmag; + else temp_dbl = 0; + scale3(temp_dbl, history); } // update history // tangential force // see e.g. eq. 18 of Thornton et al, Pow. Tech. 2013, v223,p30-46 temp_dbl = k * dt; - MathExtra::scale3(temp_dbl, vtr, temp_array); - MathExtra::sub3(history, temp_array, history); + scale3(temp_dbl, vtr, temp_array); + sub3(history, temp_array, history); } // tangential forces = history + tangential velocity damping - temp_double = -gamma; - MathExtra::scale3(temp_double, vtr, fs); + temp_dbl = -gamma_norm; + scale3(temp_dbl, vtr, fs); // rescale frictional displacements and forces if needed - fsmag = MathExtra::len3(fs); - if (fsmag > Fscrit) { - shrmag = MathExtra::len3(history); + magfs = len3(fs); + if (magfs > Fscrit) { + shrmag = len3(history); if (shrmag != 0.0) { - temp_double = Fscrit / magfs; - MathExtra::scale3(temp_double, fs, history); - MathExtra::scale3(gamma, vtr, temp_array); - MathExtra::add3(history, temp_array, history); - temp_double = Fscrit / magfs; - MathExtra::scale3(temp_double, fs); + temp_dbl = Fscrit / magfs; + scale3(temp_dbl, fs, history); + scale3(gamma_norm, vtr, temp_array); + add3(history, temp_array, history); + temp_dbl = Fscrit / magfs; + scale3(temp_dbl, fs); } else { - MathExtra::zero3(fs); + zero3(fs); } } } /* ---------------------------------------------------------------------- */ -double ContactModel::tangential_mindlin(double *history) +void ContactModel::tangential_mindlin(double *history) { - double k_scaled, gamma_scaled, fsmag, rsht, shrmag, prjmag, temp_dbl; + double k_scaled, gamma_scaled, magfs, rsht, shrmag, prjmag, temp_dbl; double temp_array[3]; int frame_update = 0; gamma_scaled = gamma_tang * damp_normal_prefactor; - k_scaled = k_tange * a; + k_scaled = k_tang * a; if (mindlin_rescale) { // on unloading, rescale the shear displacements/force if (a < history[3]) { - temp_double = a / history[3]; - MathExtra::scale3(temp_double, history); + temp_dbl = a / history[3]; + scale3(temp_dbl, history); } } // rotate and update displacements / force. // see e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235 if (history_update) { - rsht = MathExtra::dot(history, nx); + rsht = dot3(history, nx); if (mindlin_force) - frameupdate = fabs(rsht) > EPSILON * Fscrit; + frame_update = fabs(rsht) > EPSILON * Fscrit; else - frameupdate = fabs(rsht) * k_scaled > EPSILON * Fscrit; + frame_update = fabs(rsht) * k_scaled > EPSILON * Fscrit; - if (frameupdate) { - shrmag = MathExtra::len3(history); + if (frame_update) { + shrmag = len3(history); // projection - MathExtra::scale3(rsht, nx, history); + scale3(rsht, nx, history); // also rescale to preserve magnitude - prjmag = MathExtra::len3(history); - if (prjmag > 0) temp_double = shrmag / prjmag; - else temp_double = 0; - MathExtra::scale3(temp_double, history); + prjmag = len3(history); + if (prjmag > 0) temp_dbl = shrmag / prjmag; + else temp_dbl = 0; + scale3(temp_dbl, history); } // update history @@ -468,51 +478,51 @@ double ContactModel::tangential_mindlin(double *history) // tangential force // see e.g. eq. 18 of Thornton et al, Pow. Tech. 2013, v223,p30-46 temp_dbl = -k_scaled * dt; - MathExtra::scale3(temp_dbl, vtr, temp_array); + scale3(temp_dbl, vtr, temp_array); } else { - MathExtra::scale3(dt, vtr, temp_array); + scale3(dt, vtr, temp_array); } - MathExtra::add3(history, temp_array, history); + add3(history, temp_array, history); if (mindlin_rescale) history[3] = a; } // tangential forces = history + tangential velocity damping - temp_double = -gamma_scaled; - MathExtra::scale3(temp_double, vtr, fs); + temp_dbl = -gamma_scaled; + scale3(temp_dbl, vtr, fs); if (! mindlin_force) { - MathExtra::scale3(k_scaled, history, temp_array); - MathExtra::add3(fs, temp_array, fs); + scale3(k_scaled, history, temp_array); + add3(fs, temp_array, fs); } // rescale frictional displacements and forces if needed - fsmag = MathExtra::len3(fs); - if (fsmag > Fscrit) { - shrmag = MathExtra::len3(history); + magfs = len3(fs); + if (magfs > Fscrit) { + shrmag = len3(history); if (shrmag != 0.0) { - temp_double = Fscrit / magfs; - MathExtra::scale3(temp_double, fs, history); - MathExtra::scale3(gamma, vtr, temp_array); - MathExtra::add3(history, temp_array, history); + temp_dbl = Fscrit / magfs; + scale3(temp_dbl, fs, history); + scale3(gamma_tang, vtr, temp_array); + add3(history, temp_array, history); if (! mindlin_force) { - temp_double = -1.0 / k; - MathExtra::scale3(temp_double, history); + temp_dbl = -1.0 / k_tang; + scale3(temp_dbl, history); } - temp_double = Fscrit / magfs; - MathExtra::scale3(temp_double, fs); + temp_dbl = Fscrit / magfs; + scale3(temp_dbl, fs); } else { - MathExtra::zero3(fs); + zero3(fs); } } } /* ---------------------------------------------------------------------- */ -double ContactModel::rolling(double *history) +void ContactModel::rolling(double *history) { int rhist0, rhist1, rhist2, frameupdate; - double rolldotn, rollmag, prjmag, magfr, hist_temp[3], temp_array[3]; + double rolldotn, rollmag, prjmag, magfr, hist_temp[3], temp_dbl, temp_array[3]; rhist0 = roll_history_index; rhist1 = rhist0 + 1; @@ -524,46 +534,47 @@ double ContactModel::rolling(double *history) hist_temp[0] = history[rhist0]; hist_temp[1] = history[rhist1]; hist_temp[2] = history[rhist2]; - rolldotn = MathExtra::dot3(hist_temp, nx); + rolldotn = dot3(hist_temp, nx); frameupdate = fabs(rolldotn)*k_roll > EPSILON*Frcrit; if (frameupdate) { // rotate into tangential plane - rollmag = MathExtra::len3(temp); + rollmag = len3(hist_temp); // projection - temp_double = -rolldotn; - MathExtra::scale3(temp_double, nx, temp); + temp_dbl = -rolldotn; + scale3(temp_dbl, nx, temp_array); + sub3(hist_temp, temp_array, hist_temp); // also rescale to preserve magnitude - prjmag = MathExtra::len3(hist_temp); - if (prjmag > 0) temp_double = rollmag / prjmag; - else temp_double = 0; - MathExtra::scale3(temp_double, hist_temp); + prjmag = len3(hist_temp); + if (prjmag > 0) temp_dbl = rollmag / prjmag; + else temp_dbl = 0; + scale3(temp_dbl, hist_temp); } - MathExtra::scale3(dt, vrl, temp_array); - MathExtra::add3(hist_temp, temp_array, hist_temp) + scale3(dt, vrl, temp_array); + add3(hist_temp, temp_array, hist_temp); } - MathExtra::scaleadd3(k_roll, hist_tepm, gamma_roll, vrl, fr); - MathExtra::negate3(fr); + scaleadd3(k_roll, hist_temp, gamma_roll, vrl, fr); + negate3(fr); // rescale frictional displacements and forces if needed - magfr = MathExtra::len3(fr); + magfr = len3(fr); if (magfr > Frcrit) { - rollmag = MathExtra::len3(temp); + rollmag = len3(hist_temp); if (rollmag != 0.0) { - temp_double = -Frcrit / (magfr * k_roll); - MathExtra::scale3(temp_double, fr, temp_array); - MathExtra::add3(hist_temp, temp_array, hist_temp); + temp_dbl = -Frcrit / (magfr * k_roll); + scale3(temp_dbl, fr, temp_array); + add3(hist_temp, temp_array, hist_temp); - temp_double = -gamma_roll/k_roll; - MathExtra::scale3(temp_double, vrl, temp_array); - MathExtra::add3(hist_temp, temp_array, hist_temp) + temp_dbl = -gamma_roll/k_roll; + scale3(temp_dbl, vrl, temp_array); + add3(hist_temp, temp_array, hist_temp); - temp_double = Frcrit / magfr; - MathExtra::scale3(temp_double, fr); + temp_dbl = Frcrit / magfr; + scale3(temp_dbl, fr); } else { - MathExtra::zero3(fr); + zero3(fr); } } @@ -574,27 +585,27 @@ double ContactModel::rolling(double *history) /* ---------------------------------------------------------------------- */ -double ContactModel::twisting_marshall(double *history) +void ContactModel::twisting_marshall(double *history) { // Overwrite twist coefficients with derived values - k_twist = 0.5 * k_tangential * a * a; // eq 32 of Marshall paper - gamma_twist = 0.5 * gamma_tangential * a * a; - mu_twist = TWOTHIRDS * a * mu_tangential; + k_twist = 0.5 * k_tang * a * a; // eq 32 of Marshall paper + gamma_twist = 0.5 * gamma_tang * a * a; + mu_twist = TWOTHIRDS * a * mu_tang; twisting_SDS(history); } /* ---------------------------------------------------------------------- */ -double ContactModel::twisting_SDS(double *history) +void ContactModel::twisting_SDS(double *history) { double signtwist, Mtcrit; - if (historyupdate) { + if (history_update) { history[twist_history_index] += magtwist * dt; } - magtortwist = -k_twist * history[twist_history_index] - damp_twist*magtwist; // M_t torque (eq 30) + magtortwist = -k_twist * history[twist_history_index] - gamma_twist*magtwist; // M_t torque (eq 30) signtwist = (magtwist > 0) - (magtwist < 0); Mtcrit = mu_twist * Fncrit; // critical torque (eq 44) if (fabs(magtortwist) > Mtcrit) { @@ -616,9 +627,9 @@ double ContactModel::pulloff_distance(double radi, double radj) Reff_tmp = radi * radj / (radi + radj); if (Reff_tmp <= 0) return 0; - Ecaled = E * THREEQUARTERS; + Ecaled = k_norm * THREEQUARTERS; a_tmp = cbrt(9 * MY_PI * cohesion * Reff_tmp * Reff_tmp / (4 * Ecaled)); return a_tmp * a_tmp / Reff_tmp - 2 * sqrt(MY_PI * cohesion * a_tmp / Ecaled); } -}} +} diff --git a/src/GRANULAR/contact.h b/src/GRANULAR/contact.h index a5b5d23463..167f88e9de 100644 --- a/src/GRANULAR/contact.h +++ b/src/GRANULAR/contact.h @@ -11,14 +11,29 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#ifndef LMP_CONTACT -_H -#define LMP_CONTACT -_H +#ifndef LMP_CONTACT_H +#define LMP_CONTACT_H -namespace LAMMPS_NS { namespace Contact { + enum {HOOKE, HERTZ, HERTZ_MATERIAL, DMT, JKR}; + enum {VELOCITY, MASS_VELOCITY, VISCOELASTIC, TSUJI}; + enum {TANGENTIAL_NOHISTORY, TANGENTIAL_HISTORY, + TANGENTIAL_MINDLIN, TANGENTIAL_MINDLIN_RESCALE, + TANGENTIAL_MINDLIN_FORCE, TANGENTIAL_MINDLIN_RESCALE_FORCE}; + enum {TWIST_NONE, TWIST_SDS, TWIST_MARSHALL}; + enum {ROLL_NONE, ROLL_SDS}; + + #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 (4.0/3.0) // 4/3 + #define ONETHIRD (1.0/3.0) // 1/3 + #define THREEQUARTERS 0.75 // 3/4 + + #define EPSILON 1e-10 + class ContactModel { public: ContactModel(); @@ -26,6 +41,7 @@ namespace Contact { bool check_contact(); void prep_contact(); void calculate_forces(double *, double *, double *, double *); + double calculate_heat(); double pulloff_distance(double, double); int normal_model, damping_model, tangential_model; @@ -37,23 +53,25 @@ namespace Contact { double k_tang, gamma_tang, mu_tang; // tangential_coeffs - wutang? double k_roll, gamma_roll, mu_roll; // roll_coeffs double k_twist, gamma_twist, mu_twist; // twist_coeffs + double conductivity; - double radi, radj, meff, dt; - double xi[3], xj[3], vi[3], vj[3], omegai[3], omegaj[3]; + double radi, radj, meff, dt, Ti, Tj; + double *xi, *xj, *vi, *vj, *omegai, *omegaj; int history_update, roll_history_index, twist_history_index; + double fs[3], fr[3], ft[3], magtortwist; + private: - double a, knfac, Fncrit, Fscrit, Frcrit, damp_normal_prefactor; - double fs[3], fr[3], ft[3]; - double dx[3], nx[3], r, rsq, rinv, Reff, radsum, delta; - double vr[3], vn[3], vnnr, vt[3], wr[3], vtr[3], vrel; - double magtwist, magtortwist; + double a, knfac, Fntot, Fncrit, Fscrit, Frcrit, damp_normal_prefactor; + double dx[3], nx[3], r, rsq, rinv, Reff, radsum, delta, dR; + double vr[3], vn[3], vnnr, vt[3], wr[3], vtr[3], vrl[3], relrot[3], vrel; + double magtwist; bool touch; int prep_flag, check_flag; int mindlin_rescale, mindlin_force; - void touch_JKR(int); + bool touch_JKR(int); double normal_JKR(); double normal_DMT(); double normal_Hertz(); @@ -69,5 +87,5 @@ namespace Contact { }; } // namespace Contact -} // namespace LAMMPS_NS + #endif diff --git a/src/GRANULAR/fix_temp_integrate.cpp b/src/GRANULAR/fix_temp_integrate.cpp index bb351d7db4..f81442766b 100644 --- a/src/GRANULAR/fix_temp_integrate.cpp +++ b/src/GRANULAR/fix_temp_integrate.cpp @@ -36,27 +36,21 @@ FixTempIntegrate::FixTempIntegrate(LAMMPS *lmp, int narg, char **arg) : cp_style = NONE; int ntypes = atom->ntypes; - int iarg = 3; - while (iarg < narg) { - if (strcmp(arg[iarg+1],"constant") == 0) { - if (iarg+2 >= narg) error->all(FLERR,"Illegal fix command"); - cp_style = CONSTANT; - cp = utils::numeric(FLERR,arg[iarg+2],false,lmp); - if (cp < 0.0) error->all(FLERR,"Illegal fix command"); - iarg += 2; - } else if (strcmp(arg[iarg+1],"type") == 0) { - if (iarg+1+ntypes >= narg) error->all(FLERR,"Illegal fix command"); - cp_style = TYPE; - memory->create(cp_type,ntypes+1,"fix/temp/integrate:cp_type"); - for (int i = 1; i <= ntypes; i++) { - cp_type[i] = utils::numeric(FLERR,arg[iarg+1+i],false,lmp); - if (cp_type[i] < 0.0) error->all(FLERR,"Illegal fix command"); - } - iarg += 1+ntypes; - } else { - error->all(FLERR,"Illegal fix command"); + if (strcmp(arg[3],"constant") == 0) { + if (narg != 5) error->all(FLERR,"Illegal fix command"); + cp_style = CONSTANT; + cp = utils::numeric(FLERR,arg[4],false,lmp); + if (cp < 0.0) error->all(FLERR,"Illegal fix command"); + } else if (strcmp(arg[3],"type") == 0) { + if (narg != 4 + ntypes) error->all(FLERR,"Illegal fix command"); + cp_style = TYPE; + memory->create(cp_type,ntypes+1,"fix/temp/integrate:cp_type"); + for (int i = 1; i <= ntypes; i++) { + cp_type[i] = utils::numeric(FLERR,arg[3+i],false,lmp); + if (cp_type[i] < 0.0) error->all(FLERR,"Illegal fix command"); } - iarg += 1; + } else { + error->all(FLERR,"Illegal fix command"); } if (cp_style == NONE) diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index 61644df9c2..0047e24359 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -40,14 +40,13 @@ #include #include +#include using namespace LAMMPS_NS; using namespace MathConst; using namespace MathSpecial; using namespace Contact; - - /* ---------------------------------------------------------------------- */ PairGranular::PairGranular(LAMMPS *lmp) : Pair(lmp) @@ -70,10 +69,6 @@ PairGranular::PairGranular(LAMMPS *lmp) : Pair(lmp) maxrad_dynamic = nullptr; maxrad_frozen = nullptr; - models = nullptr; - - limit_damping = nullptr; - history_transfer_factors = nullptr; dt = update->dt; @@ -87,6 +82,7 @@ PairGranular::PairGranular(LAMMPS *lmp) : Pair(lmp) nondefault_history_transfer = 0; tangential_history_index = 0; roll_history_index = twist_history_index = 0; + heat_flag = 0; // create dummy fix as placeholder for FixNeighHistory // this is so final order of Modify:fix will conform to input script @@ -108,24 +104,6 @@ PairGranular::~PairGranular() if (allocated) { memory->destroy(setflag); memory->destroy(cutsq); - memory->destroy(cutoff_type); - - memory->destroy(models); - - memory->destroy(normal_coeffs); - memory->destroy(tangential_coeffs); - memory->destroy(roll_coeffs); - memory->destroy(twist_coeffs); - - memory->destroy(Emod); - memory->destroy(poiss); - - memory->destroy(normal_model); - memory->destroy(damping_model); - memory->destroy(tangential_model); - memory->destroy(roll_model); - memory->destroy(twist_model); - memory->destroy(limit_damping); delete [] onerad_dynamic; delete [] onerad_frozen; @@ -141,7 +119,7 @@ PairGranular::~PairGranular() void PairGranular::compute(int eflag, int vflag) { int i,j,ii,jj,inum,jnum,itype,jtype; - double factor_lj,mi,mj,meff; + double factor_lj,mi,mj,meff,delx,dely,delz; double forces[3], torquesi[3], torquesj[3]; int *ilist,*jlist,*numneigh,**firstneigh; @@ -184,6 +162,11 @@ void PairGranular::compute(int eflag, int vflag) int *mask = atom->mask; int nlocal = atom->nlocal; double *special_lj = force->special_lj; + double *heatflux, *temperature, dq; + if (heat_flag) { + heatflux = atom->heatflux; + temperature = atom->temperature; + } inum = list->inum; ilist = list->ilist; @@ -214,14 +197,13 @@ void PairGranular::compute(int eflag, int vflag) jtype = type[j]; // Reset model and copy initial geometric data - models[itype][jtype]->reset_contact(); + models[itype][jtype].reset_contact(); models[itype][jtype].xi = x[i]; models[itype][jtype].xj = x[j]; models[itype][jtype].radi = radius[i]; models[itype][jtype].radj = radius[j]; - - touchflag = models[itype][jtype]->check_contact(); + touchflag = models[itype][jtype].check_contact(); if (!touchflag) { // unset non-touching neighbors @@ -247,18 +229,21 @@ void PairGranular::compute(int eflag, int vflag) // Copy additional information and prepare force calculations models[itype][jtype].meff = meff; - models[itype][jtype].dt = dt; - models[itype][jtype].history_update = historyupdate; - models[itype][jtype].roll_history_index = roll_history_index; - models[itype][jtype].twist_history_index = twist_history_index; models[itype][jtype].vi = v[i]; models[itype][jtype].vj = v[j]; models[itype][jtype].omegai = omega[i]; models[itype][jtype].omegaj = omega[j]; - models[itype][jtype] -> prep_contact(); + models[itype][jtype].history_update = historyupdate; + models[itype][jtype].prep_contact(); + + if (heat_flag) { + models[itype][jtype].Ti = temperature[i]; + models[itype][jtype].Tj = temperature[j]; + } if (models[itype][jtype].normal_model == JKR) touch[jj] = 1; + // if any history is needed if (use_history) { touch[jj] = 1; @@ -266,22 +251,30 @@ void PairGranular::compute(int eflag, int vflag) } models[itype][jtype].calculate_forces(forces, torquesi, torquesj, history); + if (heat_flag) dq = models[itype][jtype].calculate_heat(); // apply forces & torques MathExtra::scale3(factor_lj, forces); MathExtra::add3(f[i], forces, f[i]); MathExtra::scale3(factor_lj, torquesi); - MathExtra::add3(torques[i], torquesi, torques[i]); + MathExtra::add3(torque[i], torquesi, torque[i]); + if (heat_flag) heatflux[i] += dq; if (force->newton_pair || j < nlocal) { MathExtra::sub3(f[j], forces, f[j]); MathExtra::scale3(factor_lj, torquesj); - MathExtra::add3(torques[j], torquesj, torques[j]); + MathExtra::add3(torque[j], torquesj, torque[j]); + if (heat_flag) heatflux[j] -= dq; } - if (evflag) ev_tally_xyz(i,j,nlocal,force->newton_pair, - 0.0,0.0,forces[0],forces[1],forces[2],dx[0],dy[1],dx[2]); + if (evflag) { + delx = x[i][0] - x[j][0]; + dely = x[i][1] - x[j][1]; + delz = x[i][2] - x[j][2]; + ev_tally_xyz(i,j,nlocal,force->newton_pair, + 0.0,0.0,forces[0],forces[1],forces[2],delx,dely,delz); + } } } } @@ -303,7 +296,9 @@ void PairGranular::allocate() memory->create(cutsq,n+1,n+1,"pair:cutsq"); memory->create(cutoff_type,n+1,n+1,"pair:cutoff_type"); - memory->create(models,n+1,n+1,"pair:models"); + for (int i = 0; i < n+1; i++) { + models.push_back(std::vector (n+1)); + } onerad_dynamic = new double[n+1]; onerad_frozen = new double[n+1]; @@ -341,6 +336,7 @@ void PairGranular::coeff(int narg, char **arg) double roll_coeffs_one[4]; double twist_coeffs_one[4]; + double conductivity_one; double cutoff_one = -1; if (narg < 2) @@ -535,6 +531,14 @@ void PairGranular::coeff(int narg, char **arg) } else if (strcmp(arg[iarg], "limit_damping") == 0) { ld_flag = 1; iarg += 1; + } else if (strcmp(arg[iarg], "heat") == 0) { + if (iarg + 1 >= narg) + error->all(FLERR, "Illegal pair_coeff command, not enough parameters"); + conductivity_one = utils::numeric(FLERR,arg[iarg+1],false,lmp); + if (conductivity_one < 0.0) + error->all(FLERR, "Illegal pair_coeff command, conductivity must be positive"); + heat_flag = 1; + iarg += 2; } else error->all(FLERR, "Illegal pair_coeff command"); } @@ -590,28 +594,29 @@ void PairGranular::coeff(int narg, char **arg) } else { models[i][j].k_tang = models[j][i].k_tang = tangential_coeffs_one[0]; } - models[i][j].gamma_tang = models[j][i]].gamma_tang = tangential_coeffs_one[1]; - models[i][j].mu_tang = models[j][i]].mu_tang = tangential_coeffs_one[2]; + models[i][j].gamma_tang = models[j][i].gamma_tang = tangential_coeffs_one[1]; + models[i][j].mu_tang = models[j][i].mu_tang = tangential_coeffs_one[2]; // Define rolling model - model[i][j].roll_model = model[j][i].roll_model = roll_model_one; + models[i][j].roll_model = models[j][i].roll_model = roll_model_one; if (roll_model_one != ROLL_NONE) { - model[i][j].k_roll = model[j][i].k_roll = roll_coeffs_one[0]; - model[i][j].gamma_roll = model[j][i].gamma_roll = roll_coeffs_one[1]; - model[i][j].mu_roll = model[j][i].mu_roll = roll_coeffs_one[2]; + models[i][j].k_roll = models[j][i].k_roll = roll_coeffs_one[0]; + models[i][j].gamma_roll = models[j][i].gamma_roll = roll_coeffs_one[1]; + models[i][j].mu_roll = models[j][i].mu_roll = roll_coeffs_one[2]; } // Define twisting model models[i][j].twist_model = models[j][i].twist_model = twist_model_one; if (twist_model_one != TWIST_NONE && twist_model_one != TWIST_MARSHALL) { - model[i][j].k_twist = model[j][i].k_twist = twist_coeffs_one[0]; - model[i][j].gamma_twist = model[j][i].gamma_twist = twist_coeffs_one[1]; - model[i][j].mu_twist = model[j][i].mu_twist = twist_coeffs_one[2]; + models[i][j].k_twist = models[j][i].k_twist = twist_coeffs_one[0]; + models[i][j].gamma_twist = models[j][i].gamma_twist = twist_coeffs_one[1]; + models[i][j].mu_twist = models[j][i].mu_twist = twist_coeffs_one[2]; } // Define extra options - model[i][j].cutoff_type = model[j][i].cutoff_type = cutoff_one; - model[i][j].limit_damping = model[j][i].limit_damping = ld_flag; + models[i][j].cutoff_type = models[j][i].cutoff_type = cutoff_one; + models[i][j].limit_damping = models[j][i].limit_damping = ld_flag; + models[i][j].conductivity = models[j][i].conductivity = conductivity_one; setflag[i][j] = 1; count++; @@ -666,8 +671,8 @@ void PairGranular::init_style() } for (int i = 1; i <= atom->ntypes; i++) for (int j = i; j <= atom->ntypes; j++) - if (model[i][j].tangential_model == TANGENTIAL_MINDLIN_RESCALE || - model[i][j].tangential_model == TANGENTIAL_MINDLIN_RESCALE_FORCE) { + if (models[i][j].tangential_model == TANGENTIAL_MINDLIN_RESCALE || + models[i][j].tangential_model == TANGENTIAL_MINDLIN_RESCALE_FORCE) { size_history += 1; roll_history_index += 1; twist_history_index += 1; @@ -785,7 +790,7 @@ double PairGranular::init_one(int i, int j) // Mix normal coefficients if (models[i][j].normal_model == HERTZ || models[i][j].normal_model == HOOKE) - models[i][j].k_norm = models[j][i][0].k_norm = + models[i][j].k_norm = models[j][i].k_norm = mix_geom(models[i][i].k_norm, models[j][j].k_norm); else models[i][j].k_norm = models[j][i].k_norm = @@ -794,7 +799,7 @@ double PairGranular::init_one(int i, int j) models[i][j].gamma_norm = models[j][i].gamma_norm = mix_geom(models[i][i].gamma_norm, models[j][j].gamma_norm); - if ((normal_model[i][j] == JKR) || (normal_model[i][j] == DMT)) + if ((models[i][j].normal_model == JKR) || (models[i][j].normal_model == DMT)) models[i][j].cohesion = models[j][i].cohesion = mix_geom(models[i][i].cohesion, models[j][j].cohesion); @@ -807,7 +812,7 @@ double PairGranular::init_one(int i, int j) mix_geom(models[i][i].mu_tang, models[j][j].mu_tang); // Mix rolling coefficients - if (models.roll_model[i][j] != ROLL_NONE) { + if (models[i][j].roll_model != ROLL_NONE) { models[i][j].k_roll = models[j][i].k_roll = mix_geom(models[i][i].k_roll, models[j][j].k_roll); models[i][j].gamma_roll = models[j][i].gamma_roll = @@ -825,8 +830,19 @@ double PairGranular::init_one(int i, int j) models[i][j].mu_twist = models[j][i].mu_twist = mix_geom(models[i][i].mu_twist, models[j][j].mu_twist); } + + models[i][j].limit_damping = models[j][i].limit_damping = + MAX(models[i][i].limit_damping, models[j][j].limit_damping); + if (heat_flag) { + models[i][j].conductivity = models[j][i].conductivity = + mix_geom(models[i][i].conductivity, models[j][j].conductivity); + } } + // Check if conductivity is undefined (negative) for any type combination + if (heat_flag && models[i][j].conductivity < 0.0) + error->all(FLERR, "To calculate heat, conductivity must be defined for all pair coeff"); + // It is possible that cut[i][j] at this point is still 0.0. // This can happen when // there is a future fix_pour after the current run. A cut[i][j] = 0.0 creates @@ -843,14 +859,14 @@ double PairGranular::init_one(int i, int j) ((maxrad_frozen[i] > 0.0) && (maxrad_dynamic[j] > 0.0))) { cutoff = maxrad_dynamic[i] + maxrad_dynamic[j]; pulloff = 0.0; - if (normal_model[i][j] == JKR) { - pulloff = models[i][j]->pulloff_distance(maxrad_dynamic[i], maxrad_dynamic[j]); + if (models[i][j].normal_model == JKR) { + pulloff = models[i][j].pulloff_distance(maxrad_dynamic[i], maxrad_dynamic[j]); cutoff += pulloff; - pulloff = models[i][j]->pulloff_distance(maxrad_frozen[i], maxrad_dynamic[j]); + pulloff = models[i][j].pulloff_distance(maxrad_frozen[i], maxrad_dynamic[j]); cutoff = MAX(cutoff, maxrad_frozen[i] + maxrad_dynamic[j] + pulloff); - pulloff = models[i][j]->pulloff_distance(maxrad_dynamic[i], maxrad_frozen[j]); + pulloff = models[i][j].pulloff_distance(maxrad_dynamic[i], maxrad_frozen[j]); cutoff = MAX(cutoff,maxrad_dynamic[i] + maxrad_frozen[j] + pulloff); } } else { @@ -872,6 +888,11 @@ double PairGranular::init_one(int i, int j) cutoff = cutoff_global; } + // Copy global options + models[i][j].dt = dt; + models[i][j].roll_history_index = roll_history_index; + models[i][j].twist_history_index = twist_history_index; + return cutoff; } @@ -982,26 +1003,114 @@ void PairGranular::reset_dt() double PairGranular::single(int i, int j, int itype, int jtype, double rsq, double /* factor_coul */, - double /* factor_lj */, double &fforce) + double factor_lj, double &fforce) { + if (factor_lj == 0) { + fforce = 0.0; + for (int m = 0; m < single_extra; m++) svector[m] = 0.0; + return 0.0; + } - fforce = Fntot*rinv; + int jnum; + int *jlist; + double *history,*allhistory; + + int nall = atom->nlocal + atom->nghost; + if ((i >= nall) || (j >= nall)) + error->all(FLERR,"Not enough atoms for pair granular single function"); + + + int touchflag; + double **x = atom->x; + double *radius = atom->radius; + + // Reset model and copy initial geometric data + models[itype][jtype].reset_contact(); + models[itype][jtype].xi = x[i]; + models[itype][jtype].xj = x[j]; + models[itype][jtype].radi = radius[i]; + models[itype][jtype].radj = radius[j]; + models[i][j].history_update = 0; // Don't update history + + touchflag = models[itype][jtype].check_contact(); + + if (!touchflag) { + fforce = 0.0; + for (int m = 0; m < single_extra; m++) svector[m] = 0.0; + return 0.0; + } + + // meff = effective mass of pair of particles + // if I or J part of rigid body, use body mass + // if I or J is frozen, meff is other particle + double mi, mj, meff; + double *rmass = atom->rmass; + int *mask = atom->mask; + + mi = rmass[i]; + mj = rmass[j]; + if (fix_rigid) { + if (mass_rigid[i] > 0.0) mi = mass_rigid[i]; + if (mass_rigid[j] > 0.0) mj = mass_rigid[j]; + } + meff = mi * mj / (mi + mj); + if (mask[i] & freeze_group_bit) meff = mj; + if (mask[j] & freeze_group_bit) meff = mi; + + // Copy additional information and prepare force calculations + double **v = atom->v; + double **omega = atom->omega; + + models[itype][jtype].meff = meff; + models[itype][jtype].vi = v[i]; + models[itype][jtype].vj = v[j]; + models[itype][jtype].omegai = omega[i]; + models[itype][jtype].omegaj = omega[j]; + models[itype][jtype].prep_contact(); + + // if any history is needed + + jnum = list->numneigh[i]; + jlist = list->firstneigh[i]; + + if (use_history) { + if ((fix_history == nullptr) || (fix_history->firstvalue == nullptr)) + error->one(FLERR,"Pair granular single computation needs history"); + allhistory = fix_history->firstvalue[i]; + for (int jj = 0; jj < jnum; jj++) { + neighprev++; + if (neighprev >= jnum) neighprev = 0; + if (jlist[neighprev] == j) break; + } + history = &allhistory[size_history*neighprev]; + } + + double forces[3], torquesi[3], torquesj[3]; + models[itype][jtype].calculate_forces(forces, torquesi, torquesj, history); + + // apply forces & torques + fforce = MathExtra::len3(forces); // set single_extra quantities - svector[0] = fs1; - svector[1] = fs2; - svector[2] = fs3; - svector[3] = fs; - svector[4] = fr1; - svector[5] = fr2; - svector[6] = fr3; - svector[7] = fr; - svector[8] = magtortwist; + double delx = x[i][0] - x[j][0]; + double dely = x[i][1] - x[j][1]; + double delz = x[i][2] - x[j][2]; + + svector[0] = models[itype][jtype].fs[0]; + svector[1] = models[itype][jtype].fs[1]; + svector[2] = models[itype][jtype].fs[2]; + svector[3] = MathExtra::len3(models[itype][jtype].fs); + svector[4] = models[itype][jtype].fr[0]; + svector[5] = models[itype][jtype].fr[1]; + svector[6] = models[itype][jtype].fr[2]; + svector[7] = MathExtra::len3(models[itype][jtype].fr); + svector[8] = models[itype][jtype].magtortwist; svector[9] = delx; svector[10] = dely; svector[11] = delz; + return 0.0; } @@ -1117,7 +1226,7 @@ double PairGranular::radii2cut(double r1, double r2) for (int i = 0; i < n; i++){ for (int j = 0; j < n; j++){ if (models[i][j].normal_model == JKR) { - temp = pulloff_distance(r1, r2); + temp = models[i][j].pulloff_distance(r1, r2); if (temp > cut) cut = temp; } } diff --git a/src/GRANULAR/pair_granular.h b/src/GRANULAR/pair_granular.h index 509981d5f9..3440d60e29 100644 --- a/src/GRANULAR/pair_granular.h +++ b/src/GRANULAR/pair_granular.h @@ -20,7 +20,9 @@ PairStyle(granular,PairGranular); #ifndef LMP_PAIR_GRANULAR_H #define LMP_PAIR_GRANULAR_H +#include "contact.h" #include "pair.h" +#include namespace LAMMPS_NS { @@ -68,9 +70,10 @@ class PairGranular : public Pair { private: int size_history; int *history_transfer_factors; + int heat_flag; // contact models - ContactModel **models; + std::vector > models; // history flags int normal_history, tangential_history, roll_history, twist_history; diff --git a/src/contact.cpp b/src/contact.cpp index b7d5d4d85d..b4b335a1aa 100644 --- a/src/contact.cpp +++ b/src/contact.cpp @@ -12,47 +12,624 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- - This class contains a series of DEM contact models - that can be defined and used to calculate forces + This class contains a series of tools for DEM contacts + Multiple models can be defined and used to calculate forces and torques based on contact geometry */ -#include #include "contact.h" +#include "math_const.h" +#include "math_extra.h" +#include "pointers.h" -namespace LAMMPS_NS { - namespace Contact{ +#include -enum {HOOKE, HERTZ, HERTZ_MATERIAL, DMT, JKR}; -enum {VELOCITY, MASS_VELOCITY, VISCOELASTIC, TSUJI}; -enum {TANGENTIAL_NOHISTORY, TANGENTIAL_HISTORY, - TANGENTIAL_MINDLIN, TANGENTIAL_MINDLIN_RESCALE, - TANGENTIAL_MINDLIN_FORCE, TANGENTIAL_MINDLIN_RESCALE_FORCE}; -enum {TWIST_NONE, TWIST_SDS, TWIST_MARSHALL}; -enum {ROLL_NONE, ROLL_SDS}; +using namespace LAMMPS_NS; +using namespace MathExtra; +using namespace MathConst; -#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 (4.0/3.0) // 4/3 -#define THREEQUARTERS 0.75 // 3/4 - -#define EPSILON 1e-10 +namespace Contact { ContactModel::ContactModel() { + k_norm = cohesion = gamma_norm = 0.0; + k_tang = gamma_tang = mu_tang = 0.0; + k_roll = gamma_roll = mu_roll = 0.0; + k_twist = gamma_twist = mu_twist = 0.0; + limit_damping = 0; + cutoff_type = 0.0; + reset_contact(); +} + +/* ---------------------------------------------------------------------- */ + +void ContactModel::reset_contact() +{ + radi = radj = 0.0; + xi = xj = vi = vj = omegai = omegaj = nullptr; + + prep_flag = check_flag = 0; + touch = false; +} + +/* ---------------------------------------------------------------------- */ + +bool ContactModel::check_contact() +{ + check_flag = 1; + + sub3(xi, xj, dx); + rsq = lensq3(dx); + radsum = radi + radj; + Reff = radi*radj/radsum; + + touch = false; + if (normal_model == JKR) touch = touch_JKR(touch); + else touch = (rsq < radsum*radsum); + + return touch; +} + +/* ---------------------------------------------------------------------- */ + +void ContactModel::prep_contact() +{ + prep_flag = 1; + + // If it hasn't already been done, test if the contact exists + if (check_flag != 1) touch = check_contact(); + if (!touch) return; + + double temp[3]; + + // Set flags + mindlin_rescale = mindlin_force = 0; + if (tangential_model == TANGENTIAL_MINDLIN_RESCALE || + tangential_model == TANGENTIAL_MINDLIN_RESCALE_FORCE) + mindlin_rescale = 1; + if (tangential_model == TANGENTIAL_MINDLIN_FORCE || + tangential_model == TANGENTIAL_MINDLIN_RESCALE_FORCE) + mindlin_force = 1; + + // Standard geometric quantities + r = sqrt(rsq); + rinv = 1.0/r; + delta = radsum - r; + dR = delta*Reff; + scale3(rinv, dx, nx); + + // relative translational velocity + sub3(vi, vj, vr); + + // normal component + vnnr = dot3(vr, nx); //v_R . n + scale3(vnnr, nx, vn); + + // tangential component + sub3(vr, vn, vt); + + // relative rotational velocity + scaleadd3(radi, omegai, radj, omegaj, wr); + + // relative tangential velocities + cross3(wr, nx, temp); + sub3(vt, temp, vtr); + vrel = len3(vtr); + + if (roll_model != ROLL_NONE || twist_model != TWIST_NONE) + sub3(omegai, omegaj, relrot); + + if (roll_model != 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) + vrl[0] = Reff * (relrot[1] * nx[2] - relrot[2] * nx[1]); + vrl[1] = Reff * (relrot[2] * nx[0] - relrot[0] * nx[2]); + vrl[2] = Reff * (relrot[0] * nx[1] - relrot[1] * nx[0]); + } + + if (twist_model != TWIST_NONE) { + // omega_T (eq 29 of Marshall) + magtwist = dot3(relrot, nx); + } +} + + +/* ---------------------------------------------------------------------- */ + +void ContactModel::calculate_forces(double *forces, double *torquesi, double *torquesj, double *history) +{ + // If it hasn't already been done, run prep calculations + if (prep_flag != 1) prep_contact(); + if (!touch) { + forces[0] = forces[1] = forces[2] = 0.0; + return; + } + + //********************************************** + // normal forces + //********************************************** + + // Also calculates: a, knfac, Fncrit (for JKR or DMT) + double Fne; + if (normal_model == JKR) { + Fne = normal_JKR(); + } else if (normal_model == DMT) { + Fne = normal_DMT(); + } else if (normal_model == HERTZ || normal_model == HERTZ_MATERIAL) { + Fne = normal_Hertz(); + } else { + Fne = normal_Hooke(); + } + + // NOTE: consider restricting Hooke to only have + // 'velocity' as an option for damping? + // Also calculates: damp_normal_prefactor + double Fdamp = normal_damping(); + + Fntot = Fne + Fdamp; + if (limit_damping && (Fntot < 0.0)) Fntot = 0.0; + + if (normal_model != JKR && normal_model != DMT) Fncrit = fabs(Fntot); + + //********************************************** + // 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 + + Fscrit = mu_tang * Fncrit; + + if (tangential_model == TANGENTIAL_MINDLIN || + tangential_model == TANGENTIAL_MINDLIN_FORCE || + tangential_model == TANGENTIAL_MINDLIN_RESCALE || + tangential_model == TANGENTIAL_MINDLIN_RESCALE_FORCE) { + tangential_mindlin(history); + } else if (tangential_model == TANGENTIAL_HISTORY) { + tangential_history(history); + } else { + tangential_no_history(); + } + + //********************************************** + // rolling force + //********************************************** + + if (roll_model != ROLL_NONE) + rolling(history); + + //**************************************** + // twisting torque, including history effects + //**************************************** + + if (twist_model == TWIST_MARSHALL) { + twisting_marshall(history); + } else if (twist_model == TWIST_SDS) { + twisting_SDS(history); + } + + //********************************************** + // sum contributions + //********************************************** + + scale3(Fntot, nx, forces); + add3(forces, fs, forces); + + cross3(nx, fs, torquesi); + copy3(torquesi, torquesj); + + double dist_to_contact = radi-0.5*delta; + scale3(dist_to_contact, torquesi); + dist_to_contact = radj-0.5*delta; + scale3(dist_to_contact, torquesj); + + double torroll[3]; + if (roll_model != ROLL_NONE) { + cross3(nx, fr, torroll); + scale3(Reff, torroll); + add3(torquesi, torroll, torquesi); + sub3(torquesj, torroll, torquesj); + } + + double tortwist[3]; + if (twist_model != TWIST_NONE) { + scale3(magtortwist, nx, tortwist); + add3(torquesi, tortwist, torquesi); + sub3(torquesj, tortwist, torquesj); + } +} + +/* ---------------------------------------------------------------------- */ + +double ContactModel::calculate_heat() +{ + double dT = Ti - Tj; + double Hc; + + //Dan is Emod the youngs modulus for all models? or do I need to scale? + Hc = 2 * conductivity * pow(FOURTHIRDS * Fntot * Reff / Emod, ONETHIRD); + + return Hc * dT; +} + +/* ---------------------------------------------------------------------- */ + +bool ContactModel::touch_JKR(int touch) +{ + double Escaled, R2, delta_pulloff, dist_pulloff; + bool touchflag; + + Escaled = k_norm * THREEQUARTERS; + if (touch) { + R2 = Reff * Reff; + a = cbrt(9.0 * MY_PI * cohesion * R2 / (4 * Escaled)); + delta_pulloff = a * a / Reff - 2 * sqrt(MY_PI * cohesion * a / Escaled); + dist_pulloff = radsum - delta_pulloff; + touchflag = (rsq < dist_pulloff * dist_pulloff); + } else { + touchflag = (rsq < radsum * radsum); + } + return touchflag; +} + +/* ---------------------------------------------------------------------- */ + +double ContactModel::normal_JKR() +{ + double Escaled, R2, dR2, t0, t1, t2, t3, t4, t5, t6; + double sqrt1, sqrt2, sqrt3, a, a2, F_pulloff, Fne; + + Escaled = k_norm * THREEQUARTERS; + + R2 = Reff * Reff; + dR2 = dR * dR; + t0 = cohesion * cohesion * R2 * R2 * Escaled; + t1 = PI27SQ*t0; + t2 = 8 * dR * dR2 * Escaled * Escaled * Escaled; + t3 = 4 * dR2 * Escaled; + + // in case sqrt(0) < 0 due to precision issues + sqrt1 = MAX(0, t0 * (t1 + 2 * t2)); + t4 = cbrt(t1 + t2 + THREEROOT3 * MY_PI * sqrt(sqrt1)); + t5 = t3 / t4 + t4 / Escaled; + sqrt2 = MAX(0, 2 * dR + t5); + t6 = sqrt(sqrt2); + sqrt3 = MAX(0, 4 * dR - t5 + SIXROOT6 * cohesion * MY_PI * R2 / (Escaled * t6)); + a = INVROOT6 * (t6 + sqrt(sqrt3)); + a2 = a * a; + Fne = Escaled * a * a2 / Reff - MY_2PI * a2 * sqrt(4 * cohesion * Escaled / (MY_PI * a)); + F_pulloff = 3 * MY_PI * cohesion * Reff; + + knfac = Escaled * a; + Fncrit = fabs(Fne + 2 * F_pulloff); + return Fne; +} + +/* ---------------------------------------------------------------------- */ + +double ContactModel::normal_DMT() +{ + a = sqrt(dR); + double Fne = a * k_norm * delta; + Fne -= 4 * MY_PI * cohesion * Reff; + double F_pulloff = 4 * MY_PI * cohesion * Reff; + + knfac = k_norm * a; + Fncrit = fabs(Fne + 2 * F_pulloff); + return Fne; +} + +/* ---------------------------------------------------------------------- */ + +double ContactModel::normal_Hertz() +{ + a = sqrt(dR); + double Fne = k_norm * delta * a; + + knfac = k_norm * a; + return Fne; +} + +/* ---------------------------------------------------------------------- */ + +double ContactModel::normal_Hooke() +{ + a = sqrt(dR); + double Fne = k_norm * delta; + + knfac = k_norm; + return Fne; +} + +/* ---------------------------------------------------------------------- */ + +double ContactModel::normal_damping() +{ + double damp_normal; + if (damping_model == VELOCITY) { + damp_normal = 1; + } else if (damping_model == MASS_VELOCITY) { + damp_normal = meff; + } else if (damping_model == VISCOELASTIC) { + damp_normal = a * meff; + } else if (damping_model == TSUJI) { + damp_normal = sqrt(meff * knfac); + } else damp_normal = 0.0; + + damp_normal_prefactor = gamma_norm * damp_normal; + return -damp_normal_prefactor * vnnr; +} + +/* ---------------------------------------------------------------------- */ + +void ContactModel::tangential_no_history() +{ + double gamma_scaled = gamma_tang * damp_normal_prefactor; + double fsmag, Ft; + + // classic pair gran/hooke (no history) + fsmag = gamma_scaled * vrel; + if (vrel != 0.0) Ft = MIN(Fscrit,fsmag) / vrel; + else Ft = 0.0; + + Ft = -Ft; + scale3(Ft, vtr, fs); +} + +/* ---------------------------------------------------------------------- */ + +void ContactModel::tangential_history(double *history) +{ + double gamma_scaled = gamma_tang * damp_normal_prefactor; + double k = k_tang; + int frame_update = 0; + double magfs, rsht, shrmag, prjmag, temp_dbl, temp_array[3]; + + // rotate and update displacements / force. + // see e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235 + if (history_update) { + rsht = dot3(history, nx); + frame_update = fabs(rsht) * k > EPSILON * Fscrit; + + if (frame_update) { + shrmag = len3(history); + // projection + scale3(rsht, nx, history); + // also rescale to preserve magnitude + prjmag = len3(history); + if (prjmag > 0) temp_dbl = shrmag / prjmag; + else temp_dbl = 0; + scale3(temp_dbl, history); + } + + // update history + // tangential force + // see e.g. eq. 18 of Thornton et al, Pow. Tech. 2013, v223,p30-46 + temp_dbl = k * dt; + scale3(temp_dbl, vtr, temp_array); + sub3(history, temp_array, history); + } + + // tangential forces = history + tangential velocity damping + temp_dbl = -gamma_norm; + scale3(temp_dbl, vtr, fs); + + // rescale frictional displacements and forces if needed + magfs = len3(fs); + if (magfs > Fscrit) { + shrmag = len3(history); + if (shrmag != 0.0) { + temp_dbl = Fscrit / magfs; + scale3(temp_dbl, fs, history); + scale3(gamma_norm, vtr, temp_array); + add3(history, temp_array, history); + temp_dbl = Fscrit / magfs; + scale3(temp_dbl, fs); + } else { + zero3(fs); + } + } +} + +/* ---------------------------------------------------------------------- */ + +void ContactModel::tangential_mindlin(double *history) +{ + double k_scaled, gamma_scaled, magfs, rsht, shrmag, prjmag, temp_dbl; + double temp_array[3]; + int frame_update = 0; + + gamma_scaled = gamma_tang * damp_normal_prefactor; + k_scaled = k_tang * a; + if (mindlin_rescale) { + // on unloading, rescale the shear displacements/force + if (a < history[3]) { + temp_dbl = a / history[3]; + scale3(temp_dbl, history); + } + } + + // rotate and update displacements / force. + // see e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235 + if (history_update) { + rsht = dot3(history, nx); + if (mindlin_force) + frame_update = fabs(rsht) > EPSILON * Fscrit; + else + frame_update = fabs(rsht) * k_scaled > EPSILON * Fscrit; + + if (frame_update) { + shrmag = len3(history); + // projection + scale3(rsht, nx, history); + // also rescale to preserve magnitude + prjmag = len3(history); + if (prjmag > 0) temp_dbl = shrmag / prjmag; + else temp_dbl = 0; + scale3(temp_dbl, history); + } + + // update history + if (mindlin_force) { + // tangential force + // see e.g. eq. 18 of Thornton et al, Pow. Tech. 2013, v223,p30-46 + temp_dbl = -k_scaled * dt; + scale3(temp_dbl, vtr, temp_array); + } else { + scale3(dt, vtr, temp_array); + } + add3(history, temp_array, history); + + if (mindlin_rescale) history[3] = a; + } + + // tangential forces = history + tangential velocity damping + temp_dbl = -gamma_scaled; + scale3(temp_dbl, vtr, fs); + + if (! mindlin_force) { + scale3(k_scaled, history, temp_array); + add3(fs, temp_array, fs); + } + + // rescale frictional displacements and forces if needed + magfs = len3(fs); + if (magfs > Fscrit) { + shrmag = len3(history); + if (shrmag != 0.0) { + temp_dbl = Fscrit / magfs; + scale3(temp_dbl, fs, history); + scale3(gamma_tang, vtr, temp_array); + add3(history, temp_array, history); + if (! mindlin_force) { + temp_dbl = -1.0 / k_tang; + scale3(temp_dbl, history); + } + temp_dbl = Fscrit / magfs; + scale3(temp_dbl, fs); + } else { + zero3(fs); + } + } +} + +/* ---------------------------------------------------------------------- */ + +void ContactModel::rolling(double *history) +{ + int rhist0, rhist1, rhist2, frameupdate; + double rolldotn, rollmag, prjmag, magfr, hist_temp[3], temp_dbl, temp_array[3]; + + rhist0 = roll_history_index; + rhist1 = rhist0 + 1; + rhist2 = rhist1 + 1; + + Frcrit = mu_roll * Fncrit; + + if (history_update) { + hist_temp[0] = history[rhist0]; + hist_temp[1] = history[rhist1]; + hist_temp[2] = history[rhist2]; + rolldotn = dot3(hist_temp, nx); + + frameupdate = fabs(rolldotn)*k_roll > EPSILON*Frcrit; + if (frameupdate) { // rotate into tangential plane + rollmag = len3(hist_temp); + // projection + temp_dbl = -rolldotn; + scale3(temp_dbl, nx, temp_array); + sub3(hist_temp, temp_array, hist_temp); + + // also rescale to preserve magnitude + prjmag = len3(hist_temp); + if (prjmag > 0) temp_dbl = rollmag / prjmag; + else temp_dbl = 0; + scale3(temp_dbl, hist_temp); + } + scale3(dt, vrl, temp_array); + add3(hist_temp, temp_array, hist_temp); + } + + scaleadd3(k_roll, hist_temp, gamma_roll, vrl, fr); + negate3(fr); + + // rescale frictional displacements and forces if needed + + magfr = len3(fr); + if (magfr > Frcrit) { + rollmag = len3(hist_temp); + if (rollmag != 0.0) { + temp_dbl = -Frcrit / (magfr * k_roll); + scale3(temp_dbl, fr, temp_array); + add3(hist_temp, temp_array, hist_temp); + + temp_dbl = -gamma_roll/k_roll; + scale3(temp_dbl, vrl, temp_array); + add3(hist_temp, temp_array, hist_temp); + + temp_dbl = Frcrit / magfr; + scale3(temp_dbl, fr); + } else { + zero3(fr); + } + } + + history[rhist0] = hist_temp[0]; + history[rhist1] = hist_temp[1]; + history[rhist2] = hist_temp[2]; +} + +/* ---------------------------------------------------------------------- */ + +void ContactModel::twisting_marshall(double *history) +{ + // Overwrite twist coefficients with derived values + k_twist = 0.5 * k_tang * a * a; // eq 32 of Marshall paper + gamma_twist = 0.5 * gamma_tang * a * a; + mu_twist = TWOTHIRDS * a * mu_tang; + + twisting_SDS(history); +} + +/* ---------------------------------------------------------------------- */ + +void ContactModel::twisting_SDS(double *history) +{ + double signtwist, Mtcrit; + + if (history_update) { + history[twist_history_index] += magtwist * dt; + } + + magtortwist = -k_twist * history[twist_history_index] - gamma_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] = (Mtcrit * signtwist - gamma_twist * magtwist) / k_twist; + magtortwist = -Mtcrit * signtwist; // eq 34 + } } /* ---------------------------------------------------------------------- - get volume-correct r basis in: basis*cbrt(vol) = q*r + compute pull-off distance (beyond contact) for a given radius and atom type + use temporary variables since this does not use a specific contact geometry ------------------------------------------------------------------------- */ -void ContactModel::() +double ContactModel::pulloff_distance(double radi, double radj) { + double Ecaled, a_tmp, Reff_tmp; + if (normal_model != JKR) return radi+radj; + + Reff_tmp = radi * radj / (radi + radj); + if (Reff_tmp <= 0) return 0; + Ecaled = k_norm * THREEQUARTERS; + a_tmp = cbrt(9 * MY_PI * cohesion * Reff_tmp * Reff_tmp / (4 * Ecaled)); + return a_tmp * a_tmp / Reff_tmp - 2 * sqrt(MY_PI * cohesion * a_tmp / Ecaled); } - -}} +} diff --git a/src/contact.h b/src/contact.h index 38cb6f0be6..167f88e9de 100644 --- a/src/contact.h +++ b/src/contact.h @@ -11,40 +11,81 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#ifndef LMP_CONTACT -_H -#define LMP_CONTACT -_H +#ifndef LMP_CONTACT_H +#define LMP_CONTACT_H -namespace LAMMPS_NS { namespace Contact { + enum {HOOKE, HERTZ, HERTZ_MATERIAL, DMT, JKR}; + enum {VELOCITY, MASS_VELOCITY, VISCOELASTIC, TSUJI}; + enum {TANGENTIAL_NOHISTORY, TANGENTIAL_HISTORY, + TANGENTIAL_MINDLIN, TANGENTIAL_MINDLIN_RESCALE, + TANGENTIAL_MINDLIN_FORCE, TANGENTIAL_MINDLIN_RESCALE_FORCE}; + enum {TWIST_NONE, TWIST_SDS, TWIST_MARSHALL}; + enum {ROLL_NONE, ROLL_SDS}; + + #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 (4.0/3.0) // 4/3 + #define ONETHIRD (1.0/3.0) // 1/3 + #define THREEQUARTERS 0.75 // 3/4 + + #define EPSILON 1e-10 + class ContactModel { public: ContactModel(); - void touch_JKR(int); - void normal_JKR(double&); - void normal_DMT(double&); - void normal_Hooke(double&); - double normal_damping(double, double); - double critical_normal(double, double); + void reset_contact(); + bool check_contact(); + void prep_contact(); + void calculate_forces(double *, double *, double *, double *); + double calculate_heat(); + double pulloff_distance(double, double); int normal_model, damping_model, tangential_model; int roll_model, twist_model; - double E, G, poisson, damp, coh; + int limit_damping; + double cutoff_type; + double Emod, poisson; // variables used in defining mixed interactions + double k_norm, gamma_norm, cohesion; // normal_coeffs + double k_tang, gamma_tang, mu_tang; // tangential_coeffs - wutang? + double k_roll, gamma_roll, mu_roll; // roll_coeffs + double k_twist, gamma_twist, mu_twist; // twist_coeffs + double conductivity; + + double radi, radj, meff, dt, Ti, Tj; + double *xi, *xj, *vi, *vj, *omegai, *omegaj; + int history_update, roll_history_index, twist_history_index; + + double fs[3], fr[3], ft[3], magtortwist; private: - double a, knfac; - ContactGeometry geom; + double a, knfac, Fntot, Fncrit, Fscrit, Frcrit, damp_normal_prefactor; + double dx[3], nx[3], r, rsq, rinv, Reff, radsum, delta, dR; + double vr[3], vn[3], vnnr, vt[3], wr[3], vtr[3], vrl[3], relrot[3], vrel; + double magtwist; + bool touch; + + int prep_flag, check_flag; + int mindlin_rescale, mindlin_force; + + bool touch_JKR(int); + double normal_JKR(); + double normal_DMT(); + double normal_Hertz(); + double normal_Hooke(); + double normal_damping(); + void tangential_no_history(); + void tangential_history(double *); + void tangential_mindlin(double *); + void rolling(double *); + void twisting_marshall(double *); + void twisting_SDS(double *); + }; - - class ContactGeometry { - public: - ContactGeometry(); - void add_data(); - double r, rinv, rsq, Reff, radsum, delta, dR; - }; } // namespace Contact -} // namespace LAMMPS_NS + #endif diff --git a/src/contact_model.cpp b/src/contact_model.cpp deleted file mode 100644 index 779631836b..0000000000 --- a/src/contact_model.cpp +++ /dev/null @@ -1,58 +0,0 @@ -// clang-format off -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - https://www.lammps.org/, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- - - This class contains a series of DEM contact models - that can be defined and used to calculate forces - and torques based on contact geometry -*/ - -#include -#include "contact_model.h" - -namespace LAMMPS_NS { - namespace Contact_Model{ - -enum {HOOKE, HERTZ, HERTZ_MATERIAL, DMT, JKR}; -enum {VELOCITY, MASS_VELOCITY, VISCOELASTIC, TSUJI}; -enum {TANGENTIAL_NOHISTORY, TANGENTIAL_HISTORY, - TANGENTIAL_MINDLIN, TANGENTIAL_MINDLIN_RESCALE, - TANGENTIAL_MINDLIN_FORCE, TANGENTIAL_MINDLIN_RESCALE_FORCE}; -enum {TWIST_NONE, TWIST_SDS, TWIST_MARSHALL}; -enum {ROLL_NONE, ROLL_SDS}; - -#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 (4.0/3.0) // 4/3 -#define THREEQUARTERS 0.75 // 3/4 - -#define EPSILON 1e-10 - -ContactModel::ContactModel() -{ - -} - -/* ---------------------------------------------------------------------- - get volume-correct r basis in: basis*cbrt(vol) = q*r -------------------------------------------------------------------------- */ - -void ContactModel::() -{ - -} - - -}} diff --git a/src/contact_model.h b/src/contact_model.h deleted file mode 100644 index 678a7111a8..0000000000 --- a/src/contact_model.h +++ /dev/null @@ -1,35 +0,0 @@ -/* -*- c++ -*- ---------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - https://www.lammps.org/, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#ifndef LMP_CONTACT_MODEL_H -#define LMP_CONTACT_MODEL_H - -namespace LAMMPS_NS { -namespace Contact_Model { - - class ContactModel { - public: - ContactModel(); - void set_strain(const double, const double); - void step_deform(const double, const double); - bool reduce(); - void get_box(double[3][3], double); - void get_rot(double[3][3]); - void get_inverse_cob(int[3][3]); - - private: - double a, knfac; - }; -} // namespace Contact_Model -} // namespace LAMMPS_NS -#endif diff --git a/src/dump_custom.cpp b/src/dump_custom.cpp index 87bd04b2a9..a74e91c889 100644 --- a/src/dump_custom.cpp +++ b/src/dump_custom.cpp @@ -41,7 +41,7 @@ enum{ID,MOL,PROC,PROCP1,TYPE,ELEMENT,MASS, XSU,YSU,ZSU,XSUTRI,YSUTRI,ZSUTRI, IX,IY,IZ, VX,VY,VZ,FX,FY,FZ, - Q,MUX,MUY,MUZ,MU,RADIUS,DIAMETER, + Q,MUX,MUY,MUZ,MU,RADIUS,DIAMETER,HEATFLUX,TEMPERATURE, OMEGAX,OMEGAY,OMEGAZ,ANGMOMX,ANGMOMY,ANGMOMZ, TQX,TQY,TQZ, COMPUTE,FIX,VARIABLE,IVEC,DVEC,IARRAY,DARRAY}; @@ -930,6 +930,18 @@ int DumpCustom::count() for (i = 0; i < nlocal; i++) dchoose[i] = 2.0*radius[i]; ptr = dchoose; nstride = 1; + } else if (thresh_array[ithresh] == HEATFLUX) { + if (!atom->heatflux_flag) + error->all(FLERR, + "Threshold for an atom property that isn't allocated"); + ptr = atom->heatflux; + nstride = 1; + } else if (thresh_array[ithresh] == TEMPERATURE) { + if (!atom->temperature_flag) + error->all(FLERR, + "Threshold for an atom property that isn't allocated"); + ptr = atom->temperature; + nstride = 1; } else if (thresh_array[ithresh] == OMEGAX) { if (!atom->omega_flag) error->all(FLERR, @@ -1380,6 +1392,16 @@ int DumpCustom::parse_fields(int narg, char **arg) error->all(FLERR,"Dumping an atom property that isn't allocated"); pack_choice[iarg] = &DumpCustom::pack_diameter; vtype[iarg] = Dump::DOUBLE; + } else if (strcmp(arg[iarg],"heatflux") == 0) { + if (!atom->heatflux_flag) + error->all(FLERR,"Dumping an atom property that isn't allocated"); + pack_choice[iarg] = &DumpCustom::pack_heatflux; + vtype[iarg] = Dump::DOUBLE; + } else if (strcmp(arg[iarg],"temperature") == 0) { + if (!atom->temperature_flag) + error->all(FLERR,"Dumping an atom property that isn't allocated"); + pack_choice[iarg] = &DumpCustom::pack_temperature; + vtype[iarg] = Dump::DOUBLE; } else if (strcmp(arg[iarg],"omegax") == 0) { if (!atom->omega_flag) error->all(FLERR,"Dumping an atom property that isn't allocated"); @@ -1845,6 +1867,8 @@ int DumpCustom::modify_param(int narg, char **arg) else if (strcmp(arg[1],"radius") == 0) thresh_array[nthresh] = RADIUS; else if (strcmp(arg[1],"diameter") == 0) thresh_array[nthresh] = DIAMETER; + else if (strcmp(arg[1],"heatflux") == 0) thresh_array[nthresh] = HEATFLUX; + else if (strcmp(arg[1],"temperature") == 0) thresh_array[nthresh] = TEMPERATURE; else if (strcmp(arg[1],"omegax") == 0) thresh_array[nthresh] = OMEGAX; else if (strcmp(arg[1],"omegay") == 0) thresh_array[nthresh] = OMEGAY; else if (strcmp(arg[1],"omegaz") == 0) thresh_array[nthresh] = OMEGAZ; @@ -2759,6 +2783,30 @@ void DumpCustom::pack_diameter(int n) /* ---------------------------------------------------------------------- */ +void DumpCustom::pack_heatflux(int n) +{ + double *heatflux = atom->heatflux; + + for (int i = 0; i < nchoose; i++) { + buf[n] = heatflux[clist[i]]; + n += size_one; + } +} + +/* ---------------------------------------------------------------------- */ + +void DumpCustom::pack_temperature(int n) +{ + double *temperature = atom->temperature; + + for (int i = 0; i < nchoose; i++) { + buf[n] = temperature[clist[i]]; + n += size_one; + } +} + +/* ---------------------------------------------------------------------- */ + void DumpCustom::pack_omegax(int n) { double **omega = atom->omega; diff --git a/src/dump_custom.h b/src/dump_custom.h index b27a9950cd..fd4fd413c3 100644 --- a/src/dump_custom.h +++ b/src/dump_custom.h @@ -193,6 +193,8 @@ class DumpCustom : public Dump { void pack_mu(int); void pack_radius(int); void pack_diameter(int); + void pack_heatflux(int); + void pack_temperature(int); void pack_omegax(int); void pack_omegay(int); diff --git a/src/fix_neigh_history.cpp b/src/fix_neigh_history.cpp index 613716c9a4..37d1b37946 100644 --- a/src/fix_neigh_history.cpp +++ b/src/fix_neigh_history.cpp @@ -353,7 +353,8 @@ void FixNeighHistory::pre_exchange_newton() // Ensure npartner is zeroed across all atoms, nall_neigh can be less than nall // when writing restarts when comm calls are made but modify->post_neighbor() isn't - for (i = 0; i < MAX(nall_neigh, atom->nall); i++) npartner[i] = 0; + int nall = atom->nlocal + atom->nghost; + for (i = 0; i < MAX(nall_neigh, nall); i++) npartner[i] = 0; tagint *tag = atom->tag; NeighList *list = pair->list; @@ -408,7 +409,7 @@ void FixNeighHistory::pre_exchange_newton() // store partner IDs and values for owned+ghost atoms // re-zero npartner to use as counter - for (i = 0; i < MAX(nall_neigh, atom->nall); i++) npartner[i] = 0; + for (i = 0; i < MAX(nall_neigh, nall); i++) npartner[i] = 0; for (ii = 0; ii < inum; ii++) { i = ilist[ii];