ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Fixing complier issues, adding dump options, initail conduction

Browse Source
This commit is contained in:
jtclemm
2022-06-18 10:11:57 -06:00
parent 468a6d9f29
commit 1faa7397d3
16 changed files with 1171 additions and 435 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
doc/src/dump.rst
View File

@ -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]
@ -107,6 +108,8 @@ Syntax
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

10
doc/src/pair_granular.rst
View File

@ -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 <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

11
src/GRANULAR/atom_vec_sphere_temp.cpp
View File

@ -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
------------------------------------------------------------------------- */

1
src/GRANULAR/atom_vec_sphere_temp.h
View File

@ -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;

383
src/GRANULAR/contact.cpp
View File

@ -17,29 +17,19 @@
and torques based on contact geometry
*/
#include <cmath>
#include "contact.h"
#include "math_const.h"
#include "math_extra.h"
#include "pointers.h"
#include <cmath>
using namespace LAMMPS_NS;
using namespace MathExtra;
using namespace MathConst;
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 THREEQUARTERS 0.75 // 3/4
#define EPSILON 1e-10
ContactModel::ContactModel()
{
k_norm = cohesion = gamma_norm = 0.0;
@ -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;
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;
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);
}
}}
}

46
src/GRANULAR/contact.h
View File

@ -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

18
src/GRANULAR/fix_temp_integrate.cpp
View File

@ -36,28 +36,22 @@ 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");
if (strcmp(arg[3],"constant") == 0) {
if (narg != 5) error->all(FLERR,"Illegal fix command");
cp_style = CONSTANT;
cp = utils::numeric(FLERR,arg[iarg+2],false,lmp);
cp = utils::numeric(FLERR,arg[4],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");
} 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[iarg+1+i],false,lmp);
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+ntypes;
} else {
error->all(FLERR,"Illegal fix command");
}
iarg += 1;
}
if (cp_style == NONE)
error->all(FLERR, "Must specify specific heat in fix temp/integrate");

249
src/GRANULAR/pair_granular.cpp
View File

@ -40,14 +40,13 @@
#include <cmath>
#include <cstring>
#include <vector>
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 <ContactModel> (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,7 +830,18 @@ 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
@ -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;
}
}

5
src/GRANULAR/pair_granular.h
View File

@ -20,7 +20,9 @@ PairStyle(granular,PairGranular);
#ifndef LMP_PAIR_GRANULAR_H
#define LMP_PAIR_GRANULAR_H
#include "contact.h"
#include "pair.h"
#include <vector>
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 <std::vector <Contact::ContactModel>> models;
// history flags
int normal_history, tangential_history, roll_history, twist_history;

627
src/contact.cpp
View File

@ -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 <cmath>
#include "contact.h"
#include "math_const.h"
#include "math_extra.h"
#include "pointers.h"
#include <cmath>
using namespace LAMMPS_NS;
using namespace MathExtra;
using namespace MathConst;
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 THREEQUARTERS 0.75 // 3/4
#define EPSILON 1e-10
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);
}
}}
}

85
src/contact.h
View File

@ -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

58
src/contact_model.cpp
View File

@ -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 <cmath>
#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::()
{
}
}}

35
src/contact_model.h
View File

@ -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

50
src/dump_custom.cpp
View File

@ -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;

2
src/dump_custom.h
View File

@ -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);

5
src/fix_neigh_history.cpp
View File

@ -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];