tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Fixing compilation errors and finishing model classes

Browse Source
This commit is contained in:
jtclemm
2022-08-01 13:12:44 -06:00
parent 8557ce7c97
commit 22de863da9
21 changed files with 903 additions and 1258 deletions
Download Patch File Download Diff File Expand all files Collapse all files

209
src/GRANULAR/contact.cpp
View File

@ -18,99 +18,140 @@
*/
#include "contact.h"
#include "pointers.h"
#include "math_const.h"
#include "math_extra.h"
#include "pointers.h"
#include "contact_sub_models.h"
#include "contact_normal_models.h"
#include "contact_tangential_models.h"
#include "contact_damping_models.h"
#include "contact_rolling_models.h"
#include "contact_twisting_models.h"
#include "contact_heat_models.h"
#include "comm.h"
#include "error.h"
#include <cmath>
using namespace LAMMPS_NS;
using namespace MathExtra;
using namespace MathConst;
using namespace LAMMPS_NS::MathConst;
using namespace LAMMPS_NS::Contact;
namespace Contact {
enum {NORMAL, TANGENTIAL, DAMPING, ROLLING, TWISTING, HEAT}
ContactModel::ContactModel()
ContactModel::ContactModel() :
Pointers(lmp)
{
limit_damping = 0;
beyond_contact = 0;
nondefault_history_transfer = 0;
cutoff_type = 0.0;
normal_model = nullptr;
tangential_model = nullptr;
damping_model = nullptr;
rolling_model = nullptr;
twisting_model = nullptr;
nmodels = 5;
reset_contact();
sub_models = {nullptr};
for (int i = 0; i < nmodels; i++) sub_models[i] = nullptr;
}
/* ---------------------------------------------------------------------- */
void ContactModel::init_model(char *model_name, int model_type)
ContactModel::~ContactModel()
{
delete[] transfer_history_factor;
}
/* ---------------------------------------------------------------------- */
void ContactModel::init_model(std::string model_name, ModelType model_type)
{
if (model_type == NORMAL) {
if (strcmp(model_name, "hooke") == 0) normal_model = new NormalHooke();
else if (strcmp(model_name, "hertz") == 0) normal_model = new NormalHertz();
else if (strcmp(model_name, "hertz/material") == 0) normal_model = new NormalHertzMaterial();
else if (strcmp(model_name, "dmt") == 0) normal_model = new NormalDMT();
else if (strcmp(model_name, "jkr") == 0) normal_model = new NormalJKR();
if (model_name == "hooke") normal_model = new NormalHooke();
else if (model_name == "hertz") normal_model = new NormalHertz();
else if (model_name == "hertz/material") normal_model = new NormalHertzMaterial();
else if (model_name == "dmt") normal_model = new NormalDMT();
else if (model_name == "jkr") normal_model = new NormalJKR();
else error->all(FLERR, "Normal model name not recognized");
sub_models[model_type] = &normal_model;
sub_models[model_type] = normal_model;
} else if (model_type == TANGENTIAL) {
if (strcmp(model_name, "linear_nohistory") == 0) tangential_model = new TangentialLinearNoHistory();
else if (strcmp(model_name, "linear_history") == 0) tangential_model = new TangentialLinearHistory();
else if (strcmp(model_name, "mindlin") == 0) tangential_model = new TangentialMindlin();
else if (strcmp(model_name, "mindlin/force") == 0) tangential_model = new TangentialMindlinForce();
else if (strcmp(model_name, "mindlin_rescale") == 0) tangential_model = new TangentialMindlinRescale();
else if (strcmp(model_name, "mindlin_rescale/force") == 0) tangential_model = new TangentialMindlinRescaleForce();
if (model_name =="linear_nohistory") tangential_model = new TangentialLinearNoHistory();
else if (model_name =="linear_history") tangential_model = new TangentialLinearHistory();
else if (model_name =="mindlin") tangential_model = new TangentialMindlin();
else if (model_name =="mindlin/force") tangential_model = new TangentialMindlinForce();
else if (model_name =="mindlin_rescale") tangential_model = new TangentialMindlinRescale();
else if (model_name =="mindlin_rescale/force") tangential_model = new TangentialMindlinRescaleForce();
else error->all(FLERR, "Tangential model name not recognized");
sub_models[model_type] = &tangential_model;
sub_models[model_type] = tangential_model;
} else if (model_type == DAMPING) {
if (strcmp(model_name, "velocity") == 0) damping_model = new DampingVelocity();
else if (strcmp(model_name, "mass_velocity") == 0) damping_model = new DampingMassVelocity();
else if (strcmp(model_name, "viscoelastic") == 0) damping_model = new DampingViscoelastic();
else if (strcmp(model_name, "tsuji") == 0) damping_model = new DampingTsuji();
if (model_name =="velocity") damping_model = new DampingVelocity();
else if (model_name =="mass_velocity") damping_model = new DampingMassVelocity();
else if (model_name =="viscoelastic") damping_model = new DampingViscoelastic();
else if (model_name =="tsuji") damping_model = new DampingTsuji();
else error->all(FLERR, "Damping model name not recognized");
sub_models[model_type] = &damping_model;
sub_models[model_type] = damping_model;
} else if (model_type == ROLLING) {
if (strcmp(model_name, "none") == 0) delete rolling_model;
else if (strcmp(model_name, "sds") == 0) rolling_model = new RollingSDS();
if (model_name =="none") delete rolling_model;
else if (model_name =="sds") rolling_model = new RollingSDS();
else error->all(FLERR, "Rolling model name not recognized");
sub_models[model_type] = &rolling_model;
sub_models[model_type] = rolling_model;
} else if (model_type == TWISTING) {
if (strcmp(model_name, "none") == 0) delete twisting_model;
else if (strcmp(model_name, "sds") == 0) twisting_model = new TwistingSDS();
else if (strcmp(model_name, "marshall") == 0) twisting_model = new TwistingMarshall();
if (model_name =="none") delete twisting_model;
else if (model_name =="sds") twisting_model = new TwistingSDS();
else if (model_name =="marshall") twisting_model = new TwistingMarshall();
else error->all(FLERR, "Twisting model name not recognized");
sub_models[model_type] = &twisting_model;
sub_models[model_type] = twisting_model;
} else if (model_type == HEAT) {
if (strcmp(model_name, "none") == 0) delete heat_model;
else if (strcmp(model_name, "area") == 0) heat_model = new HeatArea();
if (model_name =="none") delete heat_model;
else if (model_name =="area") heat_model = new HeatArea();
else error->all(FLERR, "Heat model name not recognized");
sub_models[model_type] = &heat_model;
sub_models[model_type] = heat_model;
}
sub_models[model_type]->name.assign(model_name);
sub_models[model_type]->contact = *this;
sub_models[model_type]->contact = this;
sub_models[model_type]->allocate_coeffs();
}
/* ---------------------------------------------------------------------- */
void ContactModel::init()
{
int i, j, size_cumulative;
for (i = 0; i < nmodels; i++) {
if (sub_models[i]->nondefault_history_transfer)
nondefault_history_transfer = 1;
if (sub_models[i]->beyond_contact)
beyond_contact = 1;
}
transfer_history_factor = new double(size_history);
for (i = 0; i < size_history; i++) transfer_history_factor[i] = -1.0;
if (nondefault_history_transfer) {
for (i = 0; i < size_history; i++) {
// Find which model controls this history value
size_cumulative = 0;
for (j = 0; j < nmodels; j++) {
if (size_cumulative + sub_models[j]->size_history > i) break;
size_cumulative += sub_models[j]->size_history;
}
// Check if model has nondefault transfers, if so copy its array
if (sub_models[j]->nondefault_history_transfer) {
transfer_history_factor[i] = sub_models[j]->transfer_history_factor[i - size_cumulative];
}
}
}
}
/* ---------------------------------------------------------------------- */
void ContactModel::mix_coeffs(ContactModel *c1, ContactModel *c2)
{
for (int i = 0; i < 5; i++)
sum_model[i]->mix_coeffs(c1->submodel[i], c2->submodel[i]);
for (int i = 0; i < nmodels; i++)
sub_models[i]->mix_coeffs(c1->sub_models[i], c2->sub_models[i]);
limit_damping = MAX(c1->limit_damping, c2->limit_damping);
cutoff_type = MAX(c1->cutoff_type, c2->cutoff_type);
@ -120,11 +161,17 @@ void ContactModel::mix_coeffs(ContactModel *c1, ContactModel *c2)
void ContactModel::write_restart(FILE *fp)
{
int num_char = -1;
int num_char, num_coeffs;
for (int i = 0; i < 5; i++) {
for (int i = 0; i < nmodels; i++) {
num_char = -1;
if (sub_models[i]) {
sub_models[i]->write_restart(fp);
num_char = sub_models[i]->name.length();
num_coeffs = sub_models[i]->num_coeffs;
fwrite(&num_char, sizeof(int), 1, fp);
fwrite(sub_models[i]->name.data(), sizeof(char), num_char, fp);
fwrite(&num_coeffs, sizeof(int), 1, fp);
fwrite(sub_models[i]->coeffs, sizeof(int), num_coeffs, fp);
} else {
fwrite(&num_char, sizeof(int), 1, fp);
}
@ -135,21 +182,34 @@ void ContactModel::write_restart(FILE *fp)
void ContactModel::read_restart(FILE *fp)
{
int num_char;
int num_char, num_coeff;
for (int i = 0; i < 5; i++) {
for (int i = 0; i < nmodels; i++) {
if (comm->me == 0)
utils::sfread(FLERR,&num_char,sizeof(int),1,fp,nullptr,error);
MPI_BCast(&num_char, 1, MPI_INT, 0, world);
utils::sfread(FLERR, &num_char, sizeof(int), 1, fp, nullptr, error);
MPI_Bcast(&num_char, 1, MPI_INT, 0, world);
if (num_char != -1) {
std::string model_name(num_char);
std::string model_name (num_char, ' ');
if (comm->me == 0)
utils::sfread(FLERR,const_cast<char*>(model_name.data()),sizeof(char),num_char,fp,nullptr, error);
MPI_BCast(const_cast<char*>(model_name.data()), num_char, MPI_CHAR, world);
utils::sfread(FLERR, const_cast<char*>(model_name.data()), sizeof(char),num_char, fp, nullptr, error);
MPI_Bcast(const_cast<char*>(model_name.data()), num_char, MPI_CHAR, 0, world);
init_model(const_cast<char*>(model_name.data(), i));
sub_models[i]->read_restart();
init_model(model_name, (ModelType) i);
if (comm->me == 0) {
utils::sfread(FLERR, &num_coeff, sizeof(int), 1, fp, nullptr, error);
if (num_coeff != sub_models[i]->num_coeffs)
error->one(FLERR, "Invalid contact model written to restart file");
}
MPI_Bcast(&num_coeff, 1, MPI_INT, 0, world);
if (comm->me == 0) {
utils::sfread(FLERR, sub_models[i]->coeffs, sizeof(int), num_coeff, fp, nullptr, error);
}
MPI_Bcast(sub_models[i]->coeffs, num_coeff, MPI_DOUBLE, 0, world);
sub_models[i]->coeffs_to_local();
}
}
}
@ -158,7 +218,7 @@ void ContactModel::read_restart(FILE *fp)
void ContactModel::reset_contact()
{
radi = radj = Fntot = Fncrit = magtortwist = 0.0;
radi = radj = Fntot = magtortwist = 0.0;
xi = xj = vi = vj = omegai = omegaj = nullptr;
forces = torquesi = torquesj = history = nullptr;
@ -177,10 +237,7 @@ bool ContactModel::check_contact()
radsum = radi + radj;
Reff = radi * radj / radsum;
touch = false;
if (normal_model.beyond_contact) normal_model.touch(touch);
else touch = (rsq < radsum*radsum);
touch = normal_model->touch();
return touch;
}
@ -220,7 +277,7 @@ void ContactModel::prep_contact()
sub3(vt, temp, vtr);
vrel = len3(vtr);
if (roll_model || twist_model)
if (rolling_model || twisting_model)
sub3(omegai, omegaj, relrot);
if (rolling_model) {
@ -250,20 +307,24 @@ void ContactModel::calculate_forces()
}
//**********************************************
// normal forces
// calculate forces
//**********************************************
double Fne, Fdamp;
Fne = normal_model.calculate_forces();
Fdamp = damping_model.calculate_forces();
area = normal_model->calculate_area();
Fne = normal_model->calculate_forces();
normal_model->set_knfac(); // Needed for damping
Fdamp = damping_model->calculate_forces();
normal_model->set_fncrit(); // Needed for tangential, rolling, twisting
Fntot = Fne + Fdamp;
if (limit_damping && Fntot < 0.0) Fntot = 0.0;
normal_model.set_fncrit();
tangential_model.calculate_forces();
if (rolling_model) rolling_model.calculate_forces();
if (twisting_model) twisting_model.calculate_forces();
tangential_model->calculate_forces();
if (rolling_model) rolling_model->calculate_forces();
if (twisting_model) twisting_model->calculate_forces();
//**********************************************
// sum contributions
@ -301,7 +362,7 @@ void ContactModel::calculate_forces()
double ContactModel::calculate_heat()
{
return heat_model.calculate_heat();
return heat_model->calculate_heat();
}
/* ----------------------------------------------------------------------
@ -311,7 +372,5 @@ double ContactModel::calculate_heat()
double ContactModel::pulloff_distance(double radi, double radj)
{
return normal_model.pulloff_distance(radi, radj);
}
return normal_model->pulloff_distance(radi, radj);
}

61
src/GRANULAR/contact.h
View File

@ -14,32 +14,44 @@
#ifndef LMP_CONTACT_H
#define LMP_CONTACT_H
#include "pointers.h"
#include "normal_contact_models.h"
#include "tangential_contact_models.h"
#include "damping_contact_models.h"
#include "rolling_contact_models.h"
#include "twisting_contact_models.h"
#include "heat_models.h"
using namespace LAMMPS_NS;
#include "pointers.h" // IWYU pragma: export
namespace LAMMPS_NS {
namespace Contact {
enum ModelType {
NORMAL = 0,
DAMPING = 1,
TANGENTIAL = 2,
ROLLING = 3,
TWISTING = 4,
HEAT = 5
};
#define EPSILON 1e-10
// forward declaration
class NormalModel;
class DampingModel;
class TangentialModel;
class RollingModel;
class TwistingModel;
class HeatModel;
class SubModel;
class ContactModel : protected Pointers {
public:
public:
ContactModel();
~ContactModel();
int init();
void init();
bool check_contact();
void reset_contact();
void prep_contact();
void calculate_forces();
double calculate_heat();
double pulloff_distance(double, double);
void init_model(char*, int);
void init_model(std::string, ModelType);
void mix_coeffs(ContactModel*, ContactModel*);
@ -47,21 +59,25 @@ public:
void read_restart(FILE *);
// Sub models
NormalModel normal_model;
DampingModel damping_model;
TangentialModel tangential_model;
RollingModel rolling_model;
TwistingModel twisting_model;
HeatModel heat_model;
NormalModel *normal_model;
DampingModel *damping_model;
TangentialModel *tangential_model;
RollingModel *rolling_model;
TwistingModel *twisting_model;
HeatModel *heat_model;
SubModel *sub_models[6]; // Need to resize if we add more model flavors
// Extra options
int beyond_contact, limit_damping;
int beyond_contact, limit_damping, history_update;
double cutoff_type;
// History variables
int size_history, nondefault_history_transfer;
double *transfer_history_factor;
double *history;
// Contact properties/output
double *forces, *torquesi, *torquesj;
double *history;
double radi, radj, meff, dt, Ti, Tj, area;
double Fntot, magtortwist;
@ -69,15 +85,18 @@ public:
double *xi, *xj, *vi, *vj, *omegai, *omegaj;
double fs[3], fr[3], ft[3];
private:
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;
protected:
int prep_flag, check_flag;
int nmodels;
};
} // namespace Contact
} // namespace LAMMPS_NS
#endif

40
src/GRANULAR/contact_damping_models.cpp
View File

@ -11,33 +11,22 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "damping_contact_models.h"
#include "math_const.h"
#include "contact_damping_models.h"
#include "contact_normal_models.h"
#include "contact.h"
#include <cmath>
#include "math_special.h"
using namespace MathConst;
using namespace LAMMPS_NS;
using namespace Contact;
using namespace MathSpecial;
namespace Contact{
/* ----------------------------------------------------------------------
Default damping model
------------------------------------------------------------------------- */
void DampingModel::DampingModel()
{
num_coeffs = 0;
}
/* ---------------------------------------------------------------------- */
void DampingModel::coeffs_to_local()
{
damp = contact.normal_model.damp;
damp = contact->normal_model->damp;
}
/* ----------------------------------------------------------------------
@ -46,34 +35,34 @@ void DampingModel::coeffs_to_local()
double DampingVelocity::calculate_forces()
{
return -damp * contact.vnnr;
return -damp * contact->vnnr;
}
/* ----------------------------------------------------------------------
Mass velocity damping
------------------------------------------------------------------------- */
double MassVelocity::calculate_forces()
double DampingMassVelocity::calculate_forces()
{
return -damp * contact.meff * contact.vnnr;
return -damp * contact->meff * contact->vnnr;
}
/* ----------------------------------------------------------------------
Default, viscoelastic damping
------------------------------------------------------------------------- */
double ViscoElastic::calculate_forces()
double DampingViscoelastic::calculate_forces()
{
return -damp * contact.meff * contact.area * contact.vnnr;
return -damp * contact->meff * contact->area * contact->vnnr;
}
/* ----------------------------------------------------------------------
Tsuji damping
------------------------------------------------------------------------- */
void Tsuji::coeffs_to_local()
void DampingTsuji::coeffs_to_local()
{
double tmp = contact.normal_model.damp;
double tmp = contact->normal_model->damp;
damp = 1.2728 - 4.2783 * tmp + 11.087 * square(tmp);
damp += -22.348 * cube(tmp)+ 27.467 * powint(tmp, 4);
damp += -18.022 * powint(tmp, 5) + 4.8218 * powint(tmp,6);
@ -81,8 +70,7 @@ void Tsuji::coeffs_to_local()
/* ---------------------------------------------------------------------- */
double Tsuji::calculate_forces()
double DampingTsuji::calculate_forces()
{
return -damp_scaled * sqrt(contact.meff * contact.normal_model.knfac) * contact.vnnr;
return -damp * sqrt(contact->meff * contact->normal_model->knfac) * contact->vnnr;
}

40
src/GRANULAR/contact_damping_models.h
View File

@ -11,60 +11,54 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifndef DAMPING_CONTACT_MODELS_H_
#define DAMPING_CONTACT_MODELS_H_
#ifndef CONTACT_DAMPING_MODELS_H_
#define CONTACT_DAMPING_MODELS_H_
#include "contact.h";
#include "sub_model.h"
#include "contact_sub_models.h"
namespace LAMMPS_NS {
namespace Contact {
class DampingModel:SubModel
{
class DampingModel : public SubModel {
public:
DampingModel();
virtual ~DampingModel() {};
virtual double calculate_forces() = 0;
DampingModel() {};
~DampingModel() {};
virtual void coeffs_to_local();
virtual void mix_coeffs(NormalModel*, NormalModel*);
virtual void mix_coeffs(DampingModel*, DampingModel*) {};
virtual double calculate_forces() = 0;
double damp;
};
/* ---------------------------------------------------------------------- */
class DampingVelocity:DampingModel
{
class DampingVelocity: public DampingModel {
public:
double calculate_forces();
};
/* ---------------------------------------------------------------------- */
class DampingMassVelocity:DampingModel
{
class DampingMassVelocity: public DampingModel {
public:
double calculate_forces();
};
/* ---------------------------------------------------------------------- */
class DampingViscoElastic:DampingModel
{
class DampingViscoelastic: public DampingModel {
public:
double calculate_forces();
};
/* ---------------------------------------------------------------------- */
class DampingTsuji:DampingModel
{
class DampingTsuji: public DampingModel {
public:
double calculate_forces();
void coeffs_to_local();
double calculate_forces();
};
} // namespace Contact
} // namespace LAMMPS_NS
}
#endif /*DAMPING_CONTACT_MODELS_H_ */
#endif /*CONTACT_DAMPING_MODELS_H_ */

184
src/GRANULAR/contact_heat_models.cpp
View File

@ -11,175 +11,39 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "normal_contact_models.h"
#include "math_const.h"
#include "contact_heat_models.h"
#include "contact.h"
#include <cmath>
using namespace LAMMPS_NS;
using namespace Contact;
using namespace MathConst;
/* ----------------------------------------------------------------------
Area-based heat conduction
------------------------------------------------------------------------- */
namespace Contact{
#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
// ************************
// Default behaviors where needed
// ************************
void NormalModel::set_fncrit(){
contact->Fncrit = fabs(contact->Fntot);
HeatArea::HeatArea()
{
num_coeffs = 1;
}
void NormalModel::mix_coeffs(NormalModel* imodel, NormalModel* jmodel){
for (int i = 0; i < num_coeffs; i++){
coeffs[i] = sqrt(imodel->coeffs[i]*jmodel->coeffs[i]);
}
/* ---------------------------------------------------------------------- */
void HeatArea::coeffs_to_local()
{
conductivity = coeffs[0];
}
//-----------------------------------------
/* ---------------------------------------------------------------------- */
//******************
// Hooke
//******************
void Hooke::Hooke(ContactModel &c){
contact = c;
num_coeffs = 2;
allocate_coeffs();
void HeatArea::mix_coeffs(HeatModel* imodel, HeatModel* jmodel)
{
coeffs[0] = mix_geom(imodel->coeffs[0], jmodel->coeffs[0]);
coeffs_to_local();
}
void Hooke::coeffs_to_local(){
k_norm = coeffs[0];
damp = coeffs[1];
/* ---------------------------------------------------------------------- */
double HeatArea::calculate_heat()
{
return conductivity * contact->area * (contact->Ti - contact->Tj);
}
double Hooke::calculate_forces(){
contact.area = sqrt(contact.dR);
contact.knfac = k_norm * contact.area;
return contact.knfac * contact.delta;
}
//******************
// Hertz
//******************
void Hertz::Hertz(ContactModel &c, int mat_flag){
contact = c;
material_prop_flag = mat_flag;
if (material_prop_flag){
num_coeffs = 3;
}
else{
num_coeffs = 2;
}
allocate_coeffs();
}
void Hertz::coeffs_to_local(){
if (material_prop_flag){
Emod = coeffs[0];
poiss = coeffs[1];
k_norm = 4/3*Emod;
}
else{
k_norm = coeffs[0];
}
}
double Hertz::calculate_forces(){
contact.area = sqrt(contact.dR);
contact.knfac = contact.k_norm * contact.area;
return contact.knfac * contact.delta;
}
//******************
// DMT
//******************
void DMT::DMT(ContactModel &c){
contact = c;
material_prop_flag = 1;
num_coeffs = 4;
allocate_coeffs();
}
double DMT::calculate_forces(){
contact.area = sqrt(contact.dR);
contact.knfac = contact.k_norm * contact.area;
double Fne = contact.knfac * contact.delta;
F_pulloff = 4 * MathConst::MY_PI * cohesion * contact.Reff;
Fne -= F_pulloff;
return Fne;
}
void DMT::set_fncrit(){
contact.Fncrit = fabs(contact.Fne + 2* F_pulloff);
}
//******************
// JKR
//******************
void JKR::JKR(ContactModel &c){
contact = c;
material_prop_flag = 1;
beyond_contact = 1;
num_coeffs = 4;
allocate_coeffs();
}
bool JKR::touch(int touch){
double Escaled, R2, delta_pulloff, dist_pulloff;
bool touchflag;
Escaled = k_norm * THREEQUARTERS;
if (touch) {
R2 = contact.Reff * contact.Reff;
a = cbrt(9.0 * MY_PI * cohesion * R2 / (4 * Escaled));
delta_pulloff = a * a / contact.Reff - 2 * sqrt(MY_PI * cohesion * a / Escaled);
dist_pulloff = contact.radsum - delta_pulloff;
touchflag = (contact.rsq < dist_pulloff * dist_pulloff);
} else {
touchflag = (rcontact.sq < contact.radsum * contact.radsum);
}
return touchflag;
}
double JKR::calculate_forces(){
double Escaled, R2, dR2, t0, t1, t2, t3, t4, t5, t6;
double sqrt1, sqrt2, sqrt3, 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;
return Fne;
}
void JKR::set_fncrit(){
contact.Fncrit = fabs(contact.Fne + 2 * F_pulloff);
}

83
src/GRANULAR/contact_heat_models.h
View File

@ -11,73 +11,36 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifndef NORMAL_CONTACT_MODELS_H_
#define NORMAL_CONTACT_MODELS_H_
#ifndef CONTACT_HEAT_MODELS_H_
#define CONTACT_HEAT_MODELS_H_
#include "contact.h";
#include "sub_model.h"
#include "contact_sub_models.h"
namespace LAMMPS_NS {
namespace Contact {
class NormalModel:SubModel{
public:
NormalModel(){};
virtual ~NormalModel(){};
virtual bool check_contact();
virtual void prep_contact();
virtual void set_fncrit();
virtual double calculate_forces() = 0;
virtual void coeffs_to_local();
virtual void mix_coeffs(NormalModel*, NormalModel*); //When mixing is needed
private:
int beyond_contact = 0;
int allow_limit_damping = 1;
class HeatModel : public SubModel {
public:
HeatModel() {};
~HeatModel() {};
virtual void coeffs_to_local() {};
virtual void mix_coeffs(HeatModel*, HeatModel*) {};
virtual double calculate_heat() = 0;
};
class Hooke:NormalModel{
public:
Hooke(ContactModel &c);
~Hooke(){};
/* ---------------------------------------------------------------------- */
class HeatArea: public HeatModel {
public:
HeatArea();
void coeffs_to_local();
double calculate_forces();
private:
double k_norm, damp;
void mix_coeffs(HeatModel*, HeatModel*);
double calculate_heat();
private:
double conductivity;
};
class Hertz:NormalModel{
public:
Hertz(ContactModel&, int);
~Hertz(){};
void coeffs_to_local();
double calculate_forces();
private:
double k_norm, damp, Emod, poiss;
};
class DMT:NormalModel{
public:
DMT(ContactModel &c);
~DMT(){};
void coeffs_to_local();
void coeffs_to_local(NormalModel*, NormalModel*);
double calculate_forces();
private:
double k_norm, damp, Emod, poiss, coh;
};
class JKR:NormalModel{
public:
JKR(ContactModel &c);
~JKR(){};
void coeffs_to_local();
double calculate_forces();
private:
double k_norm, damp, Emod, poiss, coh;
};
}
#endif /*NORMAL_CONTACT_MODELS_H_ */
} // namespace Contact
} // namespace LAMMPS_NS
#endif /*CONTACT_HEAT_MODELS_H_ */

184
src/GRANULAR/contact_normal_models.cpp
View File

@ -12,15 +12,13 @@
------------------------------------------------------------------------- */
#include "contact_normal_models.h"
#include "math_const.h"
#include "contact.h"
#include "math_const.h"
#include <cmath>
using namespace LAMMPS_NS;
using namespace Contact;
using namespace MathConst;
namespace Contact{
#define PI27SQ 266.47931882941264802866 // 27*PI**2
#define THREEROOT3 5.19615242270663202362 // 3*sqrt(3)
#define SIXROOT6 14.69693845669906728801 // 6*sqrt(6)
@ -40,23 +38,40 @@ NormalModel::NormalModel()
/* ---------------------------------------------------------------------- */
void NormalModel::set_fncrit()
bool NormalModel::touch()
{
Fncrit = fabs(contact.Fntot);
bool touchflag = (contact->rsq < contact->radsum * contact->radsum);
return touchflag;
}
/* ----------------------------------------------------------------------
called outside of compute(), do not assume geometry defined in contact
------------------------------------------------------------------------- */
double NormalModel::pulloff_distance(double radi, double radj)
{
return radi + radj;
}
/* ---------------------------------------------------------------------- */
void NormalModel::pulloff_distance(double radi, double radj)
double NormalModel::calculate_area()
{
return radi + radj;
return sqrt(contact->dR);
}
/* ---------------------------------------------------------------------- */
void NormalModel::set_fncrit()
{
Fncrit = fabs(contact->Fntot);
}
/* ----------------------------------------------------------------------
Hookean normal force
------------------------------------------------------------------------- */
void NormalHooke::NormalHooke()
NormalHooke::NormalHooke()
{
num_coeffs = 2;
allocate_coeffs();
@ -66,7 +81,7 @@ void NormalHooke::NormalHooke()
void NormalHooke::coeffs_to_local()
{
k_norm = coeffs[0];
k = coeffs[0];
damp = coeffs[1];
}
@ -83,17 +98,22 @@ void NormalHooke::mix_coeffs(NormalModel* imodel, NormalModel* jmodel)
double NormalHooke::calculate_forces()
{
contact.area = sqrt(contact.dR);
knfac = k_norm * contact.area;
Fne = knfac * contact.delta;
Fne = knfac * contact->delta;
return Fne;
}
/* ---------------------------------------------------------------------- */
void NormalHooke::set_knfac()
{
knfac = k * contact->area;
}
/* ----------------------------------------------------------------------
Hertzian normal force
------------------------------------------------------------------------- */
void NormalHertz::NormalHertz()
NormalHertz::NormalHertz()
{
num_coeffs = 2;
}
@ -102,7 +122,7 @@ void NormalHertz::NormalHertz()
void NormalHertz::coeffs_to_local()
{
k_norm = coeffs[0];
k = coeffs[0];
damp = coeffs[1];
}
@ -119,17 +139,23 @@ void NormalHertz::mix_coeffs(NormalModel* imodel, NormalModel* jmodel)
double NormalHertz::calculate_forces()
{
contact.area = sqrt(contact.dR);
knfac = contact.k_norm * contact.area;
Fne = knfac * contact.delta;
contact->area = sqrt(contact->dR);
Fne = knfac * contact->delta;
return Fne;
}
/* ---------------------------------------------------------------------- */
void NormalHertz::set_knfac()
{
knfac = k * contact->area;
}
/* ----------------------------------------------------------------------
Hertzian normal force with material properties
------------------------------------------------------------------------- */
void NormalHertzMaterial::NormalHertzMaterial()
NormalHertzMaterial::NormalHertzMaterial()
{
material_properties = 1;
num_coeffs = 3;
@ -142,7 +168,7 @@ void NormalHertzMaterial::coeffs_to_local()
Emod = coeffs[0];
damp = coeffs[1];
poiss = coeffs[2];
k_norm = 4 / 3 * Emod;
k = 4 / 3 * Emod;
}
/* ---------------------------------------------------------------------- */
@ -159,7 +185,7 @@ void NormalHertzMaterial::mix_coeffs(NormalModel* imodel, NormalModel* jmodel)
DMT normal force
------------------------------------------------------------------------- */
void NormalDMT::NormalDMT(ContactModel &c)
NormalDMT::NormalDMT()
{
allow_limit_damping = 0;
material_properties = 1;
@ -174,7 +200,7 @@ void NormalDMT::coeffs_to_local()
damp = coeffs[1];
poiss = coeffs[2];
cohesion = coeffs[3];
k_norm = 4 / 3 * Emod;
k = 4 / 3 * Emod;
}
/* ---------------------------------------------------------------------- */
@ -192,30 +218,35 @@ void NormalDMT::mix_coeffs(NormalModel* imodel, NormalModel* jmodel)
double NormalDMT::calculate_forces()
{
contact.area = sqrt(contact.dR);
knfac = k_norm * contact.area;
Fne = knfac * contact.delta;
F_pulloff = 4 * MathConst::MY_PI * cohesion * contact.Reff;
Fne = knfac * contact->delta;
F_pulloff = 4 * MathConst::MY_PI * cohesion * contact->Reff;
Fne -= F_pulloff;
return Fne;
}
/* ---------------------------------------------------------------------- */
void NormalDMT::set_knfac()
{
knfac = k * contact->area;
}
/* ---------------------------------------------------------------------- */
void NormalDMT::set_fncrit()
{
Fncrit = fabs(Fne + 2* F_pulloff);
Fncrit = fabs(Fne + 2 * F_pulloff);
}
/* ----------------------------------------------------------------------
JKR normal force
------------------------------------------------------------------------- */
void NormalJKR::NormalJKR(ContactModel &c)
NormalJKR::NormalJKR()
{
allow_limit_damping = 0;
material_properties = 1;
contact.beyond_contact = beyond_contact = 1;
beyond_contact = 1;
num_coeffs = 4;
}
@ -227,8 +258,8 @@ void NormalJKR::coeffs_to_local()
damp = coeffs[1];
poiss = coeffs[2];
cohesion = coeffs[3];
k_norm = 4/3*Emod;
Escaled = k_norm * THREEQUARTERS;
k = 4/3*Emod;
Escaled = k * THREEQUARTERS;
}
/* ---------------------------------------------------------------------- */
@ -244,70 +275,83 @@ void NormalJKR::mix_coeffs(NormalModel* imodel, NormalModel* jmodel)
/* ---------------------------------------------------------------------- */
bool NormalJKR::touch(int touch)
bool NormalJKR::touch()
{
double R2, delta_pulloff, dist_pulloff;
double area_at_pulloff, R2, delta_pulloff, dist_pulloff;
bool touchflag;
if (touch) {
R2 = contact.Reff * contact.Reff;
a = cbrt(9.0 * MY_PI * cohesion * R2 / (4 * Escaled));
delta_pulloff = a * a / contact.Reff - 2 * sqrt(MY_PI * cohesion * a / Escaled);
dist_pulloff = contact.radsum - delta_pulloff;
touchflag = (contact.rsq < dist_pulloff * dist_pulloff);
} else {
touchflag = (rcontact.sq < contact.radsum * contact.radsum);
}
R2 = contact->Reff * contact->Reff;
area_at_pulloff = cbrt(9.0 * MY_PI * cohesion * R2 / (4 * Escaled));
delta_pulloff = area_at_pulloff * area_at_pulloff / contact->Reff - 2 * sqrt(MY_PI * cohesion * area_at_pulloff /Escaled);
dist_pulloff = contact->radsum - delta_pulloff;
touchflag = (contact->rsq < dist_pulloff * dist_pulloff);
return touchflag;
}
/* ----------------------------------------------------------------------
called outside of compute(), do not assume geometry defined in contact
------------------------------------------------------------------------- */
double NormalJKR::pulloff_distance(double radi, double radj)
{
double area_at_pulloff, Reff_tmp;
Reff_tmp = radi * radj / (radi + radj); // May not be defined
if (Reff_tmp <= 0) return 0;
area_at_pulloff = cbrt(9 * MY_PI * cohesion * Reff_tmp * Reff_tmp / (4 * Escaled));
return area_at_pulloff * area_at_pulloff / Reff_tmp - 2 * sqrt(MY_PI * cohesion * area_at_pulloff / Escaled);
}
/* ---------------------------------------------------------------------- */
double NormalJKR::calculate_forces()
double NormalJKR::calculate_area()
{
double R2, dR2, t0, t1, t2, t3, t4, t5, t6;
double sqrt1, sqrt2, sqrt3, a2;
double sqrt1, sqrt2, sqrt3;
R2 = Reff * Reff;
dR2 = dR * dR;
R2 = contact->Reff * contact->Reff;
dR2 = contact->dR * contact->dR;
t0 = cohesion * cohesion * R2 * R2 * Escaled;
t1 = PI27SQ*t0;
t2 = 8 * dR * dR2 * Escaled * Escaled * Escaled;
t1 = PI27SQ * t0;
t2 = 8 * contact->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);
sqrt2 = MAX(0, 2 * contact->dR + t5);
t6 = sqrt(sqrt2);
sqrt3 = MAX(0, 4 * dR - t5 + SIXROOT6 * cohesion * MY_PI * R2 / (Escaled * t6));
contact.area = INVROOT6 * (t6 + sqrt(sqrt3));
a2 = contact.area * contact.area;
Fne = Escaled * contact.area * a2 / Reff - MY_2PI * a2 * sqrt(4 * cohesion * Escaled / (MY_PI * contact.area));
F_pulloff = 3 * MY_PI * cohesion * Reff;
sqrt3 = MAX(0, 4 * contact->dR - t5 + SIXROOT6 * cohesion * MY_PI * R2 / (Escaled * t6));
return INVROOT6 * (t6 + sqrt(sqrt3));
}
/* ---------------------------------------------------------------------- */
double NormalJKR::calculate_forces()
{
double a2;
a2 = contact->area * contact->area;
Fne = Escaled * contact->area * a2 / contact->Reff - MY_2PI * a2 * sqrt(4 * cohesion * Escaled / (MY_PI * contact->area));
F_pulloff = 3 * MY_PI * cohesion * contact->Reff;
knfac = Escaled * contact.area;
return Fne;
}
/* ---------------------------------------------------------------------- */
void NormalJKR::set_knfac()
{
knfac = Escaled * contact->area;
}
/* ---------------------------------------------------------------------- */
void NormalJKR::set_fncrit()
{
Fncrit = fabs(Fne + 2 * F_pulloff);
}
/* ---------------------------------------------------------------------- */
void NormalJKR::pulloff_distance(double radi, double radj)
{
double a_tmp, Reff_tmp;
Reff_tmp = radi * radj / (radi + radj);
if (Reff_tmp <= 0) return 0;
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);
}

68
src/GRANULAR/contact_normal_models.h
View File

@ -11,102 +11,108 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifndef NORMAL_CONTACT_MODELS_H_
#define NORMAL_CONTACT_MODELS_H_
#ifndef CONTACT_NORMAL_MODELS_H_
#define CONTACT_NORMAL_MODELS_H_
#include "contact.h";
#include "sub_model.h"
#include "contact_sub_models.h"
namespace LAMMPS_NS {
namespace Contact {
class NormalModel:SubModel{
class NormalModel : public SubModel {
public:
NormalModel();
virtual ~NormalModel() {};
virtual bool check_contact();
virtual void set_fncrit();
~NormalModel() {};
virtual void coeffs_to_local() {};
virtual void mix_coeffs(NormalModel*, NormalModel*) {};
virtual bool touch();
virtual double pulloff_distance(double, double);
virtual double calculate_area();
virtual double calculate_forces() = 0;
virtual void coeffs_to_local();
virtual void mix_coeffs(NormalModel*, NormalModel*);
void pulloff_distance(double, double);
virtual void set_knfac() = 0;
virtual void set_fncrit();
double damp; // Vestigial argument needed by damping
double Emod, poiss;
double Fncrit, Fne, knfac;
int material_properties;
};
/* ---------------------------------------------------------------------- */
class NormalHooke:NormalModel
{
class NormalHooke: public NormalModel {
public:
NormalHooke();
~NormalHooke() {};
void coeffs_to_local();
void mix_coeffs(NormalModel*, NormalModel*);
double calculate_forces();
double Emod, poiss;
void set_knfac();
private:
double k;
};
/* ---------------------------------------------------------------------- */
class NormalHertz:NormalModel
{
class NormalHertz: public NormalModel {
public:
NormalHertz();
~NormalHertz() {};
void coeffs_to_local();
void mix_coeffs(NormalModel*, NormalModel*);
double calculate_forces();
void set_knfac();
private:
double k_norm;
double k;
};
/* ---------------------------------------------------------------------- */
class NormalHertzMaterial:NormalHertz
{
class NormalHertzMaterial: public NormalHertz {
public:
NormalHertzMaterial();
~NormalHertzMaterial() {};
void coeffs_to_local();
void mix_coeffs(NormalModel*, NormalModel*);
private:
double k_norm;
double k;
};
/* ---------------------------------------------------------------------- */
class NormalDMT:NormalModel
{
class NormalDMT: public NormalModel {
public:
NormalDMT();
~NormalDMT() {};
void coeffs_to_local();
void mix_coeffs(NormalModel*, NormalModel*);
void set_fncrit();
double calculate_forces();
void set_knfac();
void set_fncrit();
private:
double k_norm, cohesion;
double k, cohesion;
double F_pulloff;
};
/* ---------------------------------------------------------------------- */
class NormalJKR:NormalModel
{
class NormalJKR: public NormalModel {
public:
NormalJKR();
~NormalJKR() {};
void coeffs_to_local();
void mix_coeffs(NormalModel*, NormalModel*);
void set_fncrit();
bool touch();
double pulloff_distance(double, double);
double calculate_area();
double calculate_forces();
void pulloff_distance(double, double);
void set_knfac();
void set_fncrit();
private:
double k_norm, cohesion;
double k, cohesion;
double Escaled, F_pulloff;
};
}
#endif /*NORMAL_CONTACT_MODELS_H_ */
} // namespace Contact
} // namespace LAMMPS_NS
#endif /*CONTACT_NORMAL_MODELS_H_ */

236
src/GRANULAR/contact_rolling_models.cpp
View File

@ -11,202 +11,71 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "normal_contact_models.h"
#include "contact_normal_models.h"
#include "contact_rolling_models.h"
#include "math_const.h"
#include "math_extra.h"
#include "contact.h"
#include <cmath>
using namespace LAMMPS_NS;
using namespace Contact;
using namespace MathConst;
using namespace MathExtra;
namespace Contact{
/* ----------------------------------------------------------------------
SDS rolling friction model
------------------------------------------------------------------------- */
#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
// ************************
// Default behaviors where needed
// ************************
void NormalModel::set_fncrit(){
contact->Fncrit = fabs(contact->Fntot);
RollingSDS::RollingSDS()
{
num_coeffs = 3;
size_history = 3;
}
void NormalModel::mix_coeffs(NormalModel* imodel, NormalModel* jmodel){
for (int i = 0; i < num_coeffs; i++){
coeffs[i] = sqrt(imodel->coeffs[i]*jmodel->coeffs[i]);
}
/* ---------------------------------------------------------------------- */
void RollingSDS::coeffs_to_local()
{
k = coeffs[0];
mu = coeffs[1];
gamma = coeffs[2];
}
//-----------------------------------------
/* ---------------------------------------------------------------------- */
//******************
// Hooke
//******************
void Hooke::Hooke(ContactModel &c){
contact = c;
num_coeffs = 2;
allocate_coeffs();
void RollingSDS::mix_coeffs(RollingModel* imodel, RollingModel* jmodel)
{
coeffs[0] = mix_geom(imodel->coeffs[0], jmodel->coeffs[0]);
coeffs[1] = mix_geom(imodel->coeffs[1], jmodel->coeffs[1]);
coeffs[2] = mix_geom(imodel->coeffs[2], jmodel->coeffs[2]);
coeffs_to_local();
}
void Hooke::coeffs_to_local(){
k_norm = coeffs[0];
damp = coeffs[1];
}
/* ---------------------------------------------------------------------- */
double Hooke::calculate_forces(){
contact.area = sqrt(contact.dR);
contact.knfac = k_norm * contact.area;
return contact.knfac * contact.delta;
}
//******************
// Hertz
//******************
void Hertz::Hertz(ContactModel &c, int mat_flag){
contact = c;
material_prop_flag = mat_flag;
if (material_prop_flag){
num_coeffs = 3;
}
else{
num_coeffs = 2;
}
allocate_coeffs();
}
void Hertz::coeffs_to_local(){
if (material_prop_flag){
Emod = coeffs[0];
poiss = coeffs[1];
k_norm = 4/3*Emod;
}
else{
k_norm = coeffs[0];
}
}
double Hertz::calculate_forces(){
contact.area = sqrt(contact.dR);
contact.knfac = contact.k_norm * contact.area;
return contact.knfac * contact.delta;
}
//******************
// DMT
//******************
void DMT::DMT(ContactModel &c){
contact = c;
material_prop_flag = 1;
num_coeffs = 4;
allocate_coeffs();
}
double DMT::calculate_forces(){
contact.area = sqrt(contact.dR);
contact.knfac = contact.k_norm * contact.area;
double Fne = contact.knfac * contact.delta;
F_pulloff = 4 * MathConst::MY_PI * cohesion * contact.Reff;
Fne -= F_pulloff;
return Fne;
}
void DMT::set_fncrit(){
contact.Fncrit = fabs(contact.Fne + 2* F_pulloff);
}
//******************
// JKR
//******************
void JKR::JKR(ContactModel &c){
contact = c;
material_prop_flag = 1;
beyond_contact = 1;
num_coeffs = 4;
allocate_coeffs();
}
bool JKR::touch(int touch){
double Escaled, R2, delta_pulloff, dist_pulloff;
bool touchflag;
Escaled = k_norm * THREEQUARTERS;
if (touch) {
R2 = contact.Reff * contact.Reff;
a = cbrt(9.0 * MY_PI * cohesion * R2 / (4 * Escaled));
delta_pulloff = a * a / contact.Reff - 2 * sqrt(MY_PI * cohesion * a / Escaled);
dist_pulloff = contact.radsum - delta_pulloff;
touchflag = (contact.rsq < dist_pulloff * dist_pulloff);
} else {
touchflag = (rcontact.sq < contact.radsum * contact.radsum);
}
return touchflag;
}
double JKR::calculate_forces(){
double Escaled, R2, dR2, t0, t1, t2, t3, t4, t5, t6;
double sqrt1, sqrt2, sqrt3, 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;
return Fne;
}
void JKR::set_fncrit(){
contact.Fncrit = fabs(contact.Fne + 2 * F_pulloff);
}
void ContactModel::rolling(double *history)
double RollingSDS::calculate_forces()
{
int rhist0, rhist1, rhist2, frameupdate;
double rolldotn, rollmag, prjmag, magfr, hist_temp[3], temp_dbl, temp_array[3];
double Frcrit, rolldotn, rollmag, prjmag, magfr, hist_temp[3], temp_dbl, temp_array[3];
rhist0 = roll_history_index;
rhist0 = history_index;
rhist1 = rhist0 + 1;
rhist2 = rhist1 + 1;
Frcrit = mu_roll * Fncrit;
Frcrit = mu * contact->normal_model->Fncrit;
if (history_update) {
hist_temp[0] = history[rhist0];
hist_temp[1] = history[rhist1];
hist_temp[2] = history[rhist2];
rolldotn = dot3(hist_temp, nx);
if (contact->history_update) {
hist_temp[0] = contact->history[rhist0];
hist_temp[1] = contact->history[rhist1];
hist_temp[2] = contact->history[rhist2];
rolldotn = dot3(hist_temp, contact->nx);
frameupdate = fabs(rolldotn)*k_roll > EPSILON*Frcrit;
frameupdate = fabs(rolldotn) * k > EPSILON * Frcrit;
if (frameupdate) { // rotate into tangential plane
rollmag = len3(hist_temp);
// projection
temp_dbl = -rolldotn;
scale3(temp_dbl, nx, temp_array);
scale3(temp_dbl, contact->nx, temp_array);
sub3(hist_temp, temp_array, hist_temp);
// also rescale to preserve magnitude
@ -215,35 +84,36 @@ void ContactModel::rolling(double *history)
else temp_dbl = 0;
scale3(temp_dbl, hist_temp);
}
scale3(dt, vrl, temp_array);
scale3(contact->dt, contact->vrl, temp_array);
add3(hist_temp, temp_array, hist_temp);
}
scaleadd3(k_roll, hist_temp, gamma_roll, vrl, fr);
negate3(fr);
scaleadd3(k, hist_temp, gamma, contact->vrl, contact->fr);
negate3(contact->fr);
// rescale frictional displacements and forces if needed
magfr = len3(fr);
magfr = len3(contact->fr);
if (magfr > Frcrit) {
rollmag = len3(hist_temp);
if (rollmag != 0.0) {
temp_dbl = -Frcrit / (magfr * k_roll);
scale3(temp_dbl, fr, temp_array);
temp_dbl = -Frcrit / (magfr * k);
scale3(temp_dbl, contact->fr, temp_array);
add3(hist_temp, temp_array, hist_temp);
temp_dbl = -gamma_roll/k_roll;
scale3(temp_dbl, vrl, temp_array);
temp_dbl = -gamma/k;
scale3(temp_dbl, contact->vrl, temp_array);
add3(hist_temp, temp_array, hist_temp);
temp_dbl = Frcrit / magfr;
scale3(temp_dbl, fr);
scale3(temp_dbl, contact->fr);
} else {
zero3(fr);
zero3(contact->fr);
}
}
history[rhist0] = hist_temp[0];
history[rhist1] = hist_temp[1];
history[rhist2] = hist_temp[2];
contact->history[rhist0] = hist_temp[0];
contact->history[rhist1] = hist_temp[1];
contact->history[rhist2] = hist_temp[2];
return 0;
}

79
src/GRANULAR/contact_rolling_models.h
View File

@ -11,73 +11,36 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifndef NORMAL_CONTACT_MODELS_H_
#define NORMAL_CONTACT_MODELS_H_
#ifndef CONTACT_ROLLING_MODELS_H_
#define CONTACT_ROLLING_MODELS_H_
#include "contact.h";
#include "sub_model.h"
#include "contact_sub_models.h"
namespace LAMMPS_NS {
namespace Contact {
class NormalModel:SubModel{
public:
NormalModel(){};
virtual ~NormalModel(){};
virtual bool check_contact();
virtual void prep_contact();
virtual void set_fncrit();
class RollingModel : public SubModel {
public:
RollingModel() {};
~RollingModel() {};
virtual void coeffs_to_local() {};
virtual void mix_coeffs(RollingModel*, RollingModel*) {};
virtual double calculate_forces() = 0;
virtual void coeffs_to_local();
virtual void mix_coeffs(NormalModel*, NormalModel*); //When mixing is needed
private:
int beyond_contact = 0;
int allow_limit_damping = 1;
};
class Hooke:NormalModel{
public:
Hooke(ContactModel &c);
~Hooke(){};
/* ---------------------------------------------------------------------- */
class RollingSDS: public RollingModel {
public:
RollingSDS();
void coeffs_to_local();
void mix_coeffs(RollingModel*, RollingModel*);
double calculate_forces();
private:
double k_norm, damp;
private:
double k, mu, gamma;
};
class Hertz:NormalModel{
public:
Hertz(ContactModel&, int);
~Hertz(){};
void coeffs_to_local();
double calculate_forces();
private:
double k_norm, damp, Emod, poiss;
};
class DMT:NormalModel{
public:
DMT(ContactModel &c);
~DMT(){};
void coeffs_to_local();
void coeffs_to_local(NormalModel*, NormalModel*);
double calculate_forces();
private:
double k_norm, damp, Emod, poiss, coh;
};
class JKR:NormalModel{
public:
JKR(ContactModel &c);
~JKR(){};
void coeffs_to_local();
double calculate_forces();
private:
double k_norm, damp, Emod, poiss, coh;
};
}
#endif /*NORMAL_CONTACT_MODELS_H_ */
} // namespace Contact
} // namespace LAMMPS_NS
#endif /*CONTACT_ROLLING_MODELS_H_ */

55
src/GRANULAR/contact_sub_models.cpp
View File

@ -16,23 +16,25 @@
Multiple models can be defined and used to calculate forces
and torques based on contact geometry
*/
#include "contact_sub_model.h"
#include "pointers.h"
#include "contact_sub_models.h"
#include "utils.h"
#include "error.h"
#include "comm.h"
using namespace LAMMPS_NS;
using namespace Contact;
namespace Contact{
SubModel::SubModel()
SubModel::SubModel() :
Pointers(lmp)
{
allocated = 0;
size_history = 0;
history_index = 0;
allow_limit_damping = 1;
beyond_contact = 0;
num_coeffs = 0;
nondefault_history_transfer = 0;
transfer_history_factor = nullptr;
}
/* ---------------------------------------------------------------------- */
@ -40,6 +42,7 @@ SubModel::SubModel()
SubModel::~SubModel()
{
if (allocated) delete [] coeffs;
delete [] transfer_history_factor;
}
/* ---------------------------------------------------------------------- */
@ -62,39 +65,6 @@ void SubModel::parse_coeffs(char **arg, int iarg)
coeffs_to_local();
}
/* ---------------------------------------------------------------------- */
void SubModel::write_restart(FILE *fp)
{
fwrite(&model_name.length(),sizeof(int),1,fp);
fwrite(model_name.data(),sizeof(char),model_name.length(),fp);
fwrite(&num_coeffs,sizeof(int),1,fp);
fwrite(coeffs,sizeof(int),num_coeffs,fp);
}
/* ---------------------------------------------------------------------- */
void SubModel::read_restart(FILE *fp, int num_char)
{
if (comm->me == 0) {
utils::sfread(FLERR,&num_coeffs,sizeof(int),1,fp,nullptr,error);
}
MPI_BCast(const_cast<char*>(model_name.data()), num_char, MPI_CHAR, world);
allocate_coeffs();
}
/* ---------------------------------------------------------------------- */
void SubModel::read_restart(FILE *fp)
{
int num_char;
if (me == 0) {
utils::sfread(FLERR,&num_char,sizeof(int),1,fp,nullptr,error);
}
MPI_BCast(&num_char, 1, MPI_INT, 0, world);
read_restart(fp, num_char);
}
/* ----------------------------------------------------------------------
mixing of Young's modulus (E)
------------------------------------------------------------------------- */
@ -115,7 +85,7 @@ double SubModel::mix_stiffnessG(double E1, double E2,
double pois1, double pois2)
{
double factor1 = 2 * (2 - pois1) * (1 + pois1) / E1;
double factor2 = 2 * (2 - pois2) * (1 + pois2) / E2)
double factor2 = 2 * (2 - pois2) * (1 + pois2) / E2;
return 1 / (factor1 + factor2);
}
@ -127,6 +97,3 @@ double SubModel::mix_geom(double val1, double val2)
{
return sqrt(val1 * val2);
}
}

42
src/GRANULAR/contact_sub_models.h
View File

@ -11,44 +11,48 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifndef GRANULAR_SUB_MODEL_H_
#define GRANULAR_SUB_MODEL_H_
#ifndef CONTACT_SUB_MODEL_H_
#define CONTACT_SUB_MODEL_H_
#include "contact.h"
#include "pointers.h"
using namespace LAMMPS_NS;
namespace LAMMPS_NS {
namespace Contact {
namespace Contact{
class SubModel : Pointers{
class SubModel : protected Pointers {
public:
SubModel();
virtual ~SubModel();
public:
int num_coeffs;
double *coeffs;
virtual double calculate_forces() = 0;
void read_restart();
virtual void parse_coeffs(char **, int);
void mix_coeff(SubModel*, SubModel*);
void write_restart(FILE*);
void read_restart(FILE*);
void read_restart(FILE*, int);
virtual void coeffs_to_local();
void parse_coeffs(char **, int);
virtual void mix_coeffs(SubModel*, SubModel*) {};
virtual void coeffs_to_local() {};
void allocate_coeffs();
std::string name;
private:
ContactModel &contact;
int allocated;
int size_history;
int nondefault_history_transfer;
double *transfer_history_factor;
int history_index;
int beyond_contact;
int allow_limit_damping;
ContactModel *contact;
protected:
int allocated;
double mix_stiffnessE(double, double, double, double);
double mix_stiffnessG(double, double, double, double);
double mix_geom(double, double);
};
}
} // namespace Contact
} // namespace LAMMPS_NS
#endif /* GRANULAR_SUB_MODEL_H_ */
#endif /* CONTACT_SUB_MODEL_H_ */

130
src/GRANULAR/contact_tangential_models.cpp
View File

@ -11,15 +11,18 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "tangential_contact_models.h"
#include "math_const.h"
#include "contact_damping_models.h"
#include "contact_normal_models.h"
#include "contact_tangential_models.h"
#include "contact.h"
#include "error.h"
#include "math_const.h"
#include "math_extra.h"
#include <cmath>
using namespace LAMMPS_NS;
using namespace Contact;
using namespace MathConst;
namespace Contact{
using namespace MathExtra;
/* ----------------------------------------------------------------------
Linear model with no history
@ -35,14 +38,15 @@ TangentialLinearNoHistory::TangentialLinearNoHistory()
void TangentialLinearNoHistory::coeffs_to_local()
{
k = 0.0; // No tangential stiffness with no history
xt = coeffs[0];
mu = coeffs[1];
damp = xt * contact.damping_model.damp;
damp = xt * contact->damping_model->damp;
}
/* ---------------------------------------------------------------------- */
void TangentialLinearNoHistory::mix_coeffs(NormalModel* imodel, NormalModel* jmodel)
void TangentialLinearNoHistory::mix_coeffs(TangentialModel* imodel, TangentialModel* jmodel)
{
coeffs[0] = mix_geom(imodel->coeffs[0], jmodel->coeffs[0]);
coeffs[1] = mix_geom(imodel->coeffs[1], jmodel->coeffs[1]);
@ -56,13 +60,13 @@ void TangentialLinearNoHistory::calculate_forces()
double Fscrit, fsmag, Ft;
// classic pair gran/hooke (no history)
Fscrit = mu * contact.normal_model.Fncrit
fsmag = damp * contact.vrel;
if (contact.vrel != 0.0) Ft = MIN(Fscrit, fsmag) / contact.vrel;
Fscrit = mu * contact->normal_model->Fncrit;
fsmag = damp * contact->vrel;
if (contact->vrel != 0.0) Ft = MIN(Fscrit, fsmag) / contact->vrel;
else Ft = 0.0;
Ft = -Ft;
scale3(Ft, contact.vtr, contact.fs);
scale3(Ft, contact->vtr, contact->fs);
}
/* ----------------------------------------------------------------------
@ -79,15 +83,15 @@ TangentialLinearHistory::TangentialLinearHistory()
void TangentialLinearHistory::coeffs_to_local()
{
kt = coeffs[0];
k = coeffs[0];
xt = coeffs[1];
mu = coeffs[2];
damp = xt * contact.damping_model.damp;
damp = xt * contact->damping_model->damp;
}
/* ---------------------------------------------------------------------- */
void TangentialLinearHistory::mix_coeffs(NormalModel* imodel, NormalModel* jmodel)
void TangentialLinearHistory::mix_coeffs(TangentialModel* imodel, TangentialModel* jmodel)
{
coeffs[0] = mix_geom(imodel->coeffs[0], jmodel->coeffs[0]);
coeffs[1] = mix_geom(imodel->coeffs[1], jmodel->coeffs[1]);
@ -102,20 +106,20 @@ void TangentialLinearHistory::calculate_forces()
double Fscrit, magfs, rsht, shrmag, prjmag, temp_dbl, temp_array[3];
int frame_update = 0;
Fscrit = contact.normal_model.Fncrit * mu;
Fscrit = contact->normal_model->Fncrit * mu;
double *history = & contact.history[history_index];
double *history = & contact->history[history_index];
// rotate and update displacements / force.
// see e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235
if (contact.history_update) {
rsht = dot3(history, contact.nx);
frame_update = fabs(rsht) * kt > EPSILON * Fscrit;
if (contact->history_update) {
rsht = dot3(history, contact->nx);
frame_update = fabs(rsht) * k > EPSILON * Fscrit;
if (frame_update) {
shrmag = len3(history);
// projection
scale3(rsht, contact.nx, history);
scale3(rsht, contact->nx, history);
// also rescale to preserve magnitude
prjmag = len3(history);
if (prjmag > 0) temp_dbl = shrmag / prjmag;
@ -126,28 +130,28 @@ void TangentialLinearHistory::calculate_forces()
// update history
// tangential force
// see e.g. eq. 18 of Thornton et al, Pow. Tech. 2013, v223,p30-46
temp_dbl = kt * contact.dt;
scale3(temp_dbl, contact.vtr, temp_array);
temp_dbl = k * contact->dt;
scale3(temp_dbl, contact->vtr, temp_array);
sub3(history, temp_array, history);
}
// tangential forces = history + tangential velocity damping
temp_dbl = -damp;
scale3(temp_dbl, contact.vtr, contact.fs);
scale3(temp_dbl, contact->vtr, contact->fs);
// rescale frictional displacements and forces if needed
magfs = len3(contact.fs);
magfs = len3(contact->fs);
if (magfs > Fscrit) {
shrmag = len3(history);
if (shrmag != 0.0) {
temp_dbl = Fscrit / magfs;
scale3(temp_dbl, contact.fs, history);
scale3(damp, contact.vtr, temp_array);
scale3(temp_dbl, contact->fs, history);
scale3(damp, contact->vtr, temp_array);
add3(history, temp_array, history);
temp_dbl = Fscrit / magfs;
scale3(temp_dbl, contact.fs);
scale3(temp_dbl, contact->fs);
} else {
zero3(contact.fs);
zero3(contact->fs);
}
}
}
@ -168,20 +172,22 @@ TangentialMindlin::TangentialMindlin()
void TangentialMindlin::coeffs_to_local()
{
kt = coeffs[0];
k = coeffs[0];
xt = coeffs[1];
mu = coeffs[2];
if (ke == -1)
kt = 8.0 * mix_stiffness(contact.normal_model.Emod, contact.normal_model.Emod,
contact.normal_model.poiss, contact.normal_model.poiss);
if (k == -1) {
if (!contact->normal_model->material_properties)
error->all(FLERR, "Must either specify tangential stiffness or material properties for normal model for the Mindlin tangential style");
k = 8.0 * mix_stiffnessE(contact->normal_model->Emod, contact->normal_model->Emod, contact->normal_model->poiss, contact->normal_model->poiss);
}
damp = xt * contact.damping_model.damp;
damp = xt * contact->damping_model->damp;
}
/* ---------------------------------------------------------------------- */
void TangentialMindlin::mix_coeffs(NormalModel* imodel, NormalModel* jmodel)
void TangentialMindlin::mix_coeffs(TangentialModel* imodel, TangentialModel* jmodel)
{
if (imodel->coeffs[0] == -1 || imodel->coeffs[0] == -1) coeffs[0] = -1;
else coeffs[0] = mix_geom(imodel->coeffs[0], jmodel->coeffs[0]);
@ -198,22 +204,22 @@ void TangentialMindlin::calculate_forces()
double temp_array[3];
int frame_update = 0;
double *history = & contact.history[history_index];
Fscrit = contact.normal_model.Fncrit * mu;
double *history = & contact->history[history_index];
Fscrit = contact->normal_model->Fncrit * mu;
k_scaled = k_tang * contact.area;
k_scaled = k * contact->area;
if (mindlin_rescale) {
// on unloading, rescale the shear displacements/force
if (contact.area < history[3]) {
temp_dbl = contact.area / history[3];
if (contact->area < history[3]) {
temp_dbl = contact->area / history[3];
scale3(temp_dbl, history);
}
}
// rotate and update displacements / force.
// see e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235
if (contact.history_update) {
rsht = dot3(history, contact.nx);
if (contact->history_update) {
rsht = dot3(history, contact->nx);
if (mindlin_force)
frame_update = fabs(rsht) > EPSILON * Fscrit;
else
@ -222,7 +228,7 @@ void TangentialMindlin::calculate_forces()
if (frame_update) {
shrmag = len3(history);
// projection
scale3(rsht, contact.nx, history);
scale3(rsht, contact->nx, history);
// also rescale to preserve magnitude
prjmag = len3(history);
if (prjmag > 0) temp_dbl = shrmag / prjmag;
@ -234,42 +240,42 @@ void TangentialMindlin::calculate_forces()
if (mindlin_force) {
// tangential force
// see e.g. eq. 18 of Thornton et al, Pow. Tech. 2013, v223,p30-46
temp_dbl = -k_scaled * contact.dt;
scale3(temp_dbl, contact.vtr, temp_array);
temp_dbl = -k_scaled * contact->dt;
scale3(temp_dbl, contact->vtr, temp_array);
} else {
scale3(contact.dt, contact.vtr, temp_array);
scale3(contact->dt, contact->vtr, temp_array);
}
add3(history, temp_array, history);
if (mindlin_rescale) history[3] = contact.area;
if (mindlin_rescale) history[3] = contact->area;
}
// tangential forces = history + tangential velocity damping
temp_dbl = -damp;
scale3(temp_dbl, contact.vtr, contact.fs);
scale3(temp_dbl, contact->vtr, contact->fs);
if (! mindlin_force) {
scale3(k_scaled, history, temp_array);
add3(contact.fs, temp_array, contact.fs);
add3(contact->fs, temp_array, contact->fs);
}
// rescale frictional displacements and forces if needed
magfs = len3(contact.fs);
magfs = len3(contact->fs);
if (magfs > Fscrit) {
shrmag = len3(history);
if (shrmag != 0.0) {
temp_dbl = Fscrit / magfs;
scale3(temp_dbl, contact.fs, history);
scale3(damp, contact.vtr, temp_array);
scale3(temp_dbl, contact->fs, history);
scale3(damp, contact->vtr, temp_array);
add3(history, temp_array, history);
if (! mindlin_force) {
temp_dbl = -1.0 / k_scaled;
scale3(temp_dbl, history);
}
temp_dbl = Fscrit / magfs;
scale3(temp_dbl, contact.fs);
scale3(temp_dbl, contact->fs);
} else {
zero3(contact.fs);
zero3(contact->fs);
}
}
}
@ -278,7 +284,7 @@ void TangentialMindlin::calculate_forces()
Mindlin force model
------------------------------------------------------------------------- */
void TangentialMindlinForce::TangentialMindlinForce()
TangentialMindlinForce::TangentialMindlinForce()
{
num_coeffs = 3;
size_history = 3;
@ -290,22 +296,32 @@ void TangentialMindlinForce::TangentialMindlinForce()
Mindlin rescale model
------------------------------------------------------------------------- */
void TangentialMindlinRescale::TangentialMindlinForce()
TangentialMindlinRescale::TangentialMindlinRescale()
{
num_coeffs = 3;
size_history = 4;
mindlin_force = 0;
mindlin_rescale = 1;
nondefault_history_transfer = 1;
transfer_history_factor = new double(size_history);
for (int i = 0; i < size_history; i++) transfer_history_factor[i] = -1.0;
transfer_history_factor[3] = +1;
}
/* ----------------------------------------------------------------------
Mindlin rescale force model
------------------------------------------------------------------------- */
void TangentialMindlinRescaleForce::TangentialMindlinForce()
TangentialMindlinRescaleForce::TangentialMindlinRescaleForce()
{
num_coeffs = 3;
size_history = 4;
mindlin_force = 1;
mindlin_rescale = 1;
nondefault_history_transfer = 1;
transfer_history_factor = new double(size_history);
for (int i = 0; i < size_history; i++) transfer_history_factor[i] = -1.0;
transfer_history_factor[3] = +1;
}

92
src/GRANULAR/contact_tangential_models.h
View File

@ -11,92 +11,84 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifndef TANGENTIAL_CONTACT_MODELS_H_
#define TANGENTIAL_CONTACT_MODELS_H_
#ifndef CONTACT_TANGENTIAL_MODELS_H_
#define CONTACT_TANGENTIAL_MODELS_H_
#include "contact.h";
#include "sub_model.h"
#include "contact_sub_models.h"
namespace LAMMPS_NS {
namespace Contact {
class TangentialModel:SubModel
{
public:
class TangentialModel : public SubModel {
public:
TangentialModel() {};
virtual ~TangentialModel() {};
virtual double calculate_forces() = 0;
virtual void coeffs_to_local();
virtual void mix_coeffs(TangentialModel*, TangentialModel*);
virtual void coeffs_to_local() {};
virtual void mix_coeffs(TangentialModel*, TangentialModel*) {};
virtual void calculate_forces() = 0;
int rescale_flag;
private:
int beyond_contact;
int allow_limit_damping;
double k, damp, mu; // Used by Marshall twisting model
};
/* ---------------------------------------------------------------------- */
class TangentialLinearNoHistory:TangentialModel
{
public:
TangentialLinearNoHistory()
virtual void coeffs_to_local();
virtual void mix_coeffs(TangentialModel*, TangentialModel*);
double calculate_forces();
private:
double xt, damp, mu;
class TangentialLinearNoHistory: public TangentialModel {
public:
TangentialLinearNoHistory();
void coeffs_to_local();
void mix_coeffs(TangentialModel*, TangentialModel*);
void calculate_forces();
private:
double xt;
};
/* ---------------------------------------------------------------------- */
class TangentialLinearHistory:TangentialModel
{
public:
TangentialLinearHistory()
virtual void coeffs_to_local();
virtual void mix_coeffs(TangentialModel*, TangentialModel*);
double calculate_forces();
private:
double kt, xt, damp, mu;
class TangentialLinearHistory: public TangentialModel {
public:
TangentialLinearHistory();
void coeffs_to_local();
void mix_coeffs(TangentialModel*, TangentialModel*);
void calculate_forces();
private:
double xt;
};
/* ---------------------------------------------------------------------- */
class TangentialMindlin:TangentialModel
{
public:
class TangentialMindlin: public TangentialModel {
public:
TangentialMindlin();
void coeffs_to_local();
void coeffs_to_local(TangentialModel*, TangentialModel*);
double calculate_forces();
private:
double k_norm, damp, Emod, poiss, coh;
void mix_coeffs(TangentialModel*, TangentialModel*);
void calculate_forces();
protected:
int mindlin_rescale, mindlin_force;
double xt;
};
/* ---------------------------------------------------------------------- */
class TangentialMindlinForce:TangentialMindlin
{
public:
class TangentialMindlinForce: public TangentialMindlin {
public:
TangentialMindlinForce();
};
/* ---------------------------------------------------------------------- */
class TangentialMindlinRescale:TangentialMindlin
{
public:
class TangentialMindlinRescale: public TangentialMindlin {
public:
TangentialMindlinRescale();
};
/* ---------------------------------------------------------------------- */
class TangentialMindlinRescaleForce:TangentialMindlinRescale
{
public:
TangentialMindlinForceRescale();
class TangentialMindlinRescaleForce: public TangentialMindlin {
public:
TangentialMindlinRescaleForce();
};
}
#endif /*TANGENTIAL_CONTACT_MODELS_H_ */
} // namespace Contact
} // namespace LAMMPS_NS
#endif /*CONTACT_TANGENTIAL_MODELS_H_ */

265
src/GRANULAR/contact_twisting_models.cpp
View File

@ -11,206 +11,99 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "normal_contact_models.h"
#include "math_const.h"
#include "contact_normal_models.h"
#include "contact_tangential_models.h"
#include "contact_twisting_models.h"
#include "contact.h"
#include "math_const.h"
#include <cmath>
using namespace LAMMPS_NS;
using namespace Contact;
using namespace MathConst;
namespace Contact{
/* ----------------------------------------------------------------------
Marshall twisting model
------------------------------------------------------------------------- */
#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
// ************************
// Default behaviors where needed
// ************************
void NormalModel::set_fncrit(){
contact->Fncrit = fabs(contact->Fntot);
}
void NormalModel::mix_coeffs(NormalModel* imodel, NormalModel* jmodel){
for (int i = 0; i < num_coeffs; i++){
coeffs[i] = sqrt(imodel->coeffs[i]*jmodel->coeffs[i]);
}
}
//-----------------------------------------
//******************
// Hooke
//******************
void Hooke::Hooke(ContactModel &c){
contact = c;
num_coeffs = 2;
allocate_coeffs();
}
void Hooke::coeffs_to_local(){
k_norm = coeffs[0];
damp = coeffs[1];
}
double Hooke::calculate_forces(){
contact.area = sqrt(contact.dR);
contact.knfac = k_norm * contact.area;
return contact.knfac * contact.delta;
}
//******************
// Hertz
//******************
void Hertz::Hertz(ContactModel &c, int mat_flag){
contact = c;
material_prop_flag = mat_flag;
if (material_prop_flag){
num_coeffs = 3;
}
else{
num_coeffs = 2;
}
allocate_coeffs();
}
void Hertz::coeffs_to_local(){
if (material_prop_flag){
Emod = coeffs[0];
poiss = coeffs[1];
k_norm = 4/3*Emod;
}
else{
k_norm = coeffs[0];
}
}
double Hertz::calculate_forces(){
contact.area = sqrt(contact.dR);
contact.knfac = contact.k_norm * contact.area;
return contact.knfac * contact.delta;
}
//******************
// DMT
//******************
void DMT::DMT(ContactModel &c){
contact = c;
material_prop_flag = 1;
num_coeffs = 4;
allocate_coeffs();
}
double DMT::calculate_forces(){
contact.area = sqrt(contact.dR);
contact.knfac = contact.k_norm * contact.area;
double Fne = contact.knfac * contact.delta;
F_pulloff = 4 * MathConst::MY_PI * cohesion * contact.Reff;
Fne -= F_pulloff;
return Fne;
}
void DMT::set_fncrit(){
contact.Fncrit = fabs(contact.Fne + 2* F_pulloff);
}
//******************
// JKR
//******************
void JKR::JKR(ContactModel &c){
contact = c;
material_prop_flag = 1;
beyond_contact = 1;
num_coeffs = 4;
allocate_coeffs();
}
bool JKR::touch(int touch){
double Escaled, R2, delta_pulloff, dist_pulloff;
bool touchflag;
Escaled = k_norm * THREEQUARTERS;
if (touch) {
R2 = contact.Reff * contact.Reff;
a = cbrt(9.0 * MY_PI * cohesion * R2 / (4 * Escaled));
delta_pulloff = a * a / contact.Reff - 2 * sqrt(MY_PI * cohesion * a / Escaled);
dist_pulloff = contact.radsum - delta_pulloff;
touchflag = (contact.rsq < dist_pulloff * dist_pulloff);
} else {
touchflag = (rcontact.sq < contact.radsum * contact.radsum);
}
return touchflag;
}
double JKR::calculate_forces(){
double Escaled, R2, dR2, t0, t1, t2, t3, t4, t5, t6;
double sqrt1, sqrt2, sqrt3, 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;
return Fne;
}
void JKR::set_fncrit(){
contact.Fncrit = fabs(contact.Fne + 2 * F_pulloff);
}
/* ---------------------------------------------------------------------- */
void ContactModel::twisting_marshall(double *history)
TwistingMarshall::TwistingMarshall()
{
// 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);
num_coeffs = 0;
size_history = 3;
}
/* ---------------------------------------------------------------------- */
void ContactModel::twisting_SDS(double *history)
double TwistingMarshall::calculate_forces()
{
double signtwist, Mtcrit;
if (history_update) {
history[twist_history_index] += magtwist * dt;
// Calculate twist coefficients from tangential model & contact geometry
// eq 32 of Marshall paper
double k = 0.5 * contact->tangential_model->k * contact->area * contact->area;
double damp = 0.5 * contact->tangential_model->damp * contact->area * contact->area;
double mu = TWOTHIRDS * contact->area * contact->tangential_model->mu;
if (contact->history_update) {
contact->history[history_index] += contact->magtwist * contact->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
// M_t torque (eq 30)
contact->magtortwist = -k * contact->history[history_index] - damp * contact->magtwist;
signtwist = (contact->magtwist > 0) - (contact->magtwist < 0);
Mtcrit = mu * contact->normal_model->Fncrit; // critical torque (eq 44)
if (fabs(contact->magtortwist) > Mtcrit) {
contact->history[history_index] = (Mtcrit * signtwist - damp * contact->magtwist) / k;
contact->magtortwist = -Mtcrit * signtwist; // eq 34
}
}
/* ----------------------------------------------------------------------
SDS twisting model
------------------------------------------------------------------------- */
TwistingSDS::TwistingSDS()
{
num_coeffs = 3;
size_history = 3;
}
/* ---------------------------------------------------------------------- */
void TwistingSDS::coeffs_to_local()
{
k = coeffs[0];
mu = coeffs[1];
damp = coeffs[2];
}
/* ---------------------------------------------------------------------- */
void TwistingSDS::mix_coeffs(TwistingModel* imodel, TwistingModel* jmodel)
{
coeffs[0] = mix_geom(imodel->coeffs[0], jmodel->coeffs[0]);
coeffs[1] = mix_geom(imodel->coeffs[1], jmodel->coeffs[1]);
coeffs[2] = mix_geom(imodel->coeffs[2], jmodel->coeffs[2]);
coeffs_to_local();
}
/* ---------------------------------------------------------------------- */
double TwistingSDS::calculate_forces()
{
double signtwist, Mtcrit;
if (contact->history_update) {
contact->history[history_index] += contact->magtwist * contact->dt;
}
// M_t torque (eq 30)
contact->magtortwist = -k * contact->history[history_index] - damp * contact->magtwist;
signtwist = (contact->magtwist > 0) - (contact->magtwist < 0);
Mtcrit = mu * contact->normal_model->Fncrit; // critical torque (eq 44)
if (fabs(contact->magtortwist) > Mtcrit) {
contact->history[history_index] = (Mtcrit * signtwist - damp * contact->magtwist) / k;
contact->magtortwist = -Mtcrit * signtwist; // eq 34
}
}

81
src/GRANULAR/contact_twisting_models.h
View File

@ -11,73 +11,44 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifndef NORMAL_CONTACT_MODELS_H_
#define NORMAL_CONTACT_MODELS_H_
#ifndef CONTACT_TWISTING_MODELS_H_
#define CONTACT_TWISTING_MODELS_H_
#include "contact.h";
#include "sub_model.h"
#include "contact_sub_models.h"
namespace LAMMPS_NS {
namespace Contact {
class NormalModel:SubModel{
public:
NormalModel(){};
virtual ~NormalModel(){};
virtual bool check_contact();
virtual void prep_contact();
virtual void set_fncrit();
class TwistingModel : public SubModel {
public:
TwistingModel() {};
virtual ~TwistingModel() {};
virtual void coeffs_to_local() {};
virtual void mix_coeffs(TwistingModel*, TwistingModel*) {};
virtual double calculate_forces() = 0;
virtual void coeffs_to_local();
virtual void mix_coeffs(NormalModel*, NormalModel*); //When mixing is needed
private:
int beyond_contact = 0;
int allow_limit_damping = 1;
};
class Hooke:NormalModel{
public:
Hooke(ContactModel &c);
~Hooke(){};
void coeffs_to_local();
/* ---------------------------------------------------------------------- */
class TwistingMarshall: public TwistingModel {
public:
TwistingMarshall();
double calculate_forces();
private:
double k_norm, damp;
};
class Hertz:NormalModel{
public:
Hertz(ContactModel&, int);
~Hertz(){};
/* ---------------------------------------------------------------------- */
class TwistingSDS: public TwistingModel {
public:
TwistingSDS();
void coeffs_to_local();
void mix_coeffs(TwistingModel*, TwistingModel*);
double calculate_forces();
private:
double k_norm, damp, Emod, poiss;
private:
double k, mu, damp;
};
class DMT:NormalModel{
public:
DMT(ContactModel &c);
~DMT(){};
void coeffs_to_local();
void coeffs_to_local(NormalModel*, NormalModel*);
double calculate_forces();
private:
double k_norm, damp, Emod, poiss, coh;
};
class JKR:NormalModel{
public:
JKR(ContactModel &c);
~JKR(){};
void coeffs_to_local();
double calculate_forces();
private:
double k_norm, damp, Emod, poiss, coh;
};
}
#endif /*NORMAL_CONTACT_MODELS_H_ */
} // namespace Contact
} // namespace LAMMPS_NS
#endif /*CONTACT_TWISTING_MODELS_H_ */

169
src/GRANULAR/pair_granular.cpp
View File

@ -23,14 +23,19 @@
#include "atom.h"
#include "comm.h"
#include "contact.h"
#include "contact_sub_models.h"
#include "contact_normal_models.h"
#include "contact_tangential_models.h"
#include "contact_damping_models.h"
#include "contact_rolling_models.h"
#include "contact_twisting_models.h"
#include "contact_heat_models.h"
#include "error.h"
#include "fix.h"
#include "fix_dummy.h"
#include "fix_neigh_history.h"
#include "force.h"
#include "math_const.h"
#include "math_extra.h"
#include "math_special.h"
#include "memory.h"
#include "modify.h"
#include "neigh_list.h"
@ -38,14 +43,12 @@
#include "neighbor.h"
#include "update.h"
#include <cmath>
#include <cstring>
#include <vector>
using namespace LAMMPS_NS;
using namespace MathConst;
using namespace MathSpecial;
using namespace Contact;
using namespace MathExtra;
/* ---------------------------------------------------------------------- */
@ -69,8 +72,6 @@ PairGranular::PairGranular(LAMMPS *lmp) : Pair(lmp)
maxrad_dynamic = nullptr;
maxrad_frozen = nullptr;
history_transfer_factors = nullptr;
dt = update->dt;
// set comm size needed by this Pair if used with fix rigid
@ -97,7 +98,6 @@ PairGranular::PairGranular(LAMMPS *lmp) : Pair(lmp)
PairGranular::~PairGranular()
{
delete[] svector;
delete[] history_transfer_factors;
if (!fix_history) modify->delete_fix("NEIGH_HISTORY_GRANULAR_DUMMY");
else modify->delete_fix("NEIGH_HISTORY_GRANULAR");
@ -121,7 +121,7 @@ void PairGranular::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double factor_lj,mi,mj,meff,delx,dely,delz;
double forces[3], torquesi[3], torquesj[3];
double *forces, *torquesi, *torquesj;
int *ilist,*jlist,*numneigh,**firstneigh;
int *touch,**firsttouch;
@ -259,17 +259,17 @@ void PairGranular::compute(int eflag, int vflag)
torquesj = models[itype][jtype]->torquesj;
// apply forces & torques
MathExtra::scale3(factor_lj, forces);
MathExtra::add3(f[i], forces, f[i]);
scale3(factor_lj, forces);
add3(f[i], forces, f[i]);
MathExtra::scale3(factor_lj, torquesi);
MathExtra::add3(torque[i], torquesi, torque[i]);
scale3(factor_lj, torquesi);
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(torque[j], torquesj, torque[j]);
sub3(f[j], forces, f[j]);
scale3(factor_lj, torquesj);
add3(torque[j], torquesj, torque[j]);
if (heat_flag) heatflux[j] -= dq;
}
@ -301,8 +301,7 @@ void PairGranular::allocate()
memory->create(cutsq,n+1,n+1,"pair:cutsq");
memory->create(cutoff_type,n+1,n+1,"pair:cutoff_type");
models = (ContactModel **) memory->srealloc(fix, n+1,n+1, "pair:contact_models");
memory->create(models,n+1,n+1,"pair:contact_models");
onerad_dynamic = new double[n+1];
onerad_frozen = new double[n+1];
@ -349,26 +348,26 @@ void PairGranular::coeff(int narg, char **arg)
//Parse mandatory normal and tangential specifications
int iarg = 2;
vec_models.back().init(arg[iarg], Contact::NORMAL);
ncoeff = vec_models.back().normal_model.num_coeffs
vec_models.back().init_model(std::string(arg[iarg]), NORMAL);
ncoeff = vec_models.back().normal_model->num_coeffs;
iarg += 1;
if (iarg + ncoeff >= narg)
error->all(FLERR,"Illegal pair_coeff command"
"Insufficient arguments provided for normal model.");
vec_models.back().normal_model.parse_coeffs(arg, iarg);
vec_models.back().normal_model->parse_coeffs(arg, iarg);
iarg += ncoeff;
if (strcmp(arg[iarg], "tangential") == 0) {
if (iarg + 1 >= narg)
error->all(FLERR,"Illegal pair_coeff command, must specify "
"tangential model after tangential keyword");
vec_models.back().init(arg[iarg], Contact::TANGENTIAL);
ncoeff = vec_models.back().tangential_model.num_coeffs;
vec_models.back().init_model(std::string(arg[iarg]), TANGENTIAL);
ncoeff = vec_models.back().tangential_model->num_coeffs;
iarg += 1;
if (iarg + ncoeff >= narg)
error->all(FLERR, "Illegal pair_coeff command"
"Insufficient arguments provided for tangential model.");
vec_models.back().tangential_model.parse_coeffs(arg, iarg);
vec_models.back().tangential_model->parse_coeffs(arg, iarg);
iarg += ncoeff;
} else{
error->all(FLERR, "Illegal pair_coeff command, 'tangential' keyword expected");
@ -379,49 +378,49 @@ void PairGranular::coeff(int narg, char **arg)
if (strcmp(arg[iarg], "damping") == 0) {
if (iarg + 1 >= narg)
error->all(FLERR, "Illegal pair_coeff command, not enough parameters");
vec_models.back().init(arg[iarg], Contact::DAMPING);
ncoeff = vec_models.back().damping_model.num_coeffs;
vec_models.back().init_model(std::string(arg[iarg]), DAMPING);
ncoeff = vec_models.back().damping_model->num_coeffs;
iarg += 1;
if (iarg + ncoeff >= narg)
error->all(FLERR, "Illegal pair_coeff command"
"Insufficient arguments provided for damping model.");
vec_models.back().damping_model.parse_coeffs(arg, iarg);
vec_models.back().damping_model->parse_coeffs(arg, iarg);
iarg += ncoeff;
} else if (strcmp(arg[iarg], "rolling") == 0) {
if (iarg + 1 >= narg)
error->all(FLERR, "Illegal pair_coeff command, not enough parameters");
vec_models.back().init(arg[iarg], Contact::ROLLING);
ncoeff = vec_models.back().rolling_model.num_coeffs;
vec_models.back().init_model(std::string(arg[iarg]), ROLLING);
ncoeff = vec_models.back().rolling_model->num_coeffs;
iarg += 1;
if (iarg + ncoeff >= narg)
error->all(FLERR, "Illegal pair_coeff command"
"Insufficient arguments provided for rolling model.");
vec_models.back().rolling_model.parse_coeffs(arg, iarg);
vec_models.back().rolling_model->parse_coeffs(arg, iarg);
iarg += ncoeff;
} else if (strcmp(arg[iarg], "twisting") == 0) {
if (iarg + 1 >= narg)
error->all(FLERR, "Illegal pair_coeff command, not enough parameters");
vec_models.back().init(arg[iarg], Contact::TWISTING);
ncoeff = vec_models.back().twisting_model.num_coeffs;
vec_models.back().init_model(std::string(arg[iarg]), TWISTING);
ncoeff = vec_models.back().twisting_model->num_coeffs;
iarg += 1;
if (iarg + ncoeff >= narg)
error->all(FLERR, "Illegal pair_coeff command"
"Insufficient arguments provided for twisting model.");
vec_models.back().twisting_model.parse_coeffs(arg, iarg);
vec_models.back().twisting_model->parse_coeffs(arg, iarg);
iarg += ncoeff;
} else if (strcmp(arg[iarg], "heat") == 0) {
if (iarg + 1 >= narg)
error->all(FLERR, "Illegal pair_coeff command, not enough parameters");
vec_models.back().init(arg[iarg], Contact::HEAT);
ncoeff = vec_models.back().heat_model.num_coeffs;
vec_models.back().init_model(std::string(arg[iarg]), HEAT);
ncoeff = vec_models.back().heat_model->num_coeffs;
iarg += 1;
if (iarg + ncoeff >= narg)
error->all(FLERR, "Illegal pair_coeff command"
"Insufficient arguments provided for heat model.");
vec_models.back().heat_model.parse_coeffs(arg, iarg);
vec_models.back().heat_model->parse_coeffs(arg, iarg);
iarg += ncoeff;
heat_flag = 1;
@ -440,8 +439,8 @@ void PairGranular::coeff(int narg, char **arg)
// Define default damping model if unspecified, takes no args
if (!vec_models.back().damping_model) {
vec_models.back().init("viscoelastic", Contact::DAMPING);
vec_models.back().damping_model.parse_coeffs(arg, 0);
vec_models.back().init_model("viscoelastic", DAMPING);
vec_models.back().damping_model->parse_coeffs(arg, 0);
}
if (vec_models.back().limit_damping && !vec_models.back().normal_model->allow_limit_damping)
@ -475,38 +474,41 @@ void PairGranular::init_style()
if (comm->ghost_velocity == 0)
error->all(FLERR,"Pair granular requires ghost atoms store velocity");
// allocate history
// allocate history and initialize models
int size_normal_history = 0;
int size_damping_history = 0;
int size_tangential_history = 0;
int size_rolling_history = 0;
int size_twisting_history = 0;
for (int i = 1; i <= atom->ntypes; i++) {
for (int j = i; j <= atom->ntypes; j++) {
if (models[i][j]->normal_model.history_flag ||
models[i][j]->tangential_model.history_flag ||
models[i][j]->rolling_model.history_flag ||
models[i][j]->twisting_model.history_flag) use_history = 1;
if (models[i][j]->nondefault_history_transfer) // ???
if (models[i][j]->normal_model->size_history != 0 ||
models[i][j]->damping_model->size_history != 0 ||
models[i][j]->tangential_model->size_history != 0 ||
models[i][j]->rolling_model->size_history != 0 ||
models[i][j]->twisting_model->size_history != 0) use_history = 1;
if (models[i][j]->normal_model->size_history > size_normal_history)
size_normal_history = models[i][j]->normal_model.size_history;
size_normal_history = models[i][j]->damping_model->size_history;
if (models[i][j]->damping_model->size_history > size_damping_history)
size_damping_history = models[i][j]->normal_model->size_history;
if (models[i][j]->tangential_model->size_history > size_tangential_history)
size_tangential_history = models[i][j]->tangential_model.size_history;
size_tangential_history = models[i][j]->tangential_model->size_history;
if (models[i][j]->rolling_model->size_history > size_rolling_history)
size_rolling_history = models[i][j]->rolling_model.size_history;
size_rolling_history = models[i][j]->rolling_model->size_history;
if (models[i][j]->twisting_model->size_history > size_twisting_history)
size_twisting_history = models[i][j]->twisting_model.size_history;
size_twisting_history = models[i][j]->twisting_model->size_history;
}
}
size_history = size_normal_history + size_tangential_history +
size_rolling_history + size_twisting_history;
size_history = size_normal_history + size_damping_history +
size_tangential_history + size_rolling_history + size_twisting_history;
tangential_history_index = size_normal_history;
roll_history_index = size_normal_history + size_tangential_history;
twist_history_index = size_normal_history + size_tangential_history + size_rolling_history;
damping_history_index = size_normal_history;
tangential_history_index = size_normal_history + damping_history_index;
roll_history_index = size_normal_history + damping_history_index + size_tangential_history;
twist_history_index = size_normal_history + damping_history_index + size_tangential_history + size_rolling_history;
if (use_history) neighbor->add_request(this, NeighConst::REQ_SIZE|NeighConst::REQ_HISTORY);
else neighbor->add_request(this, NeighConst::REQ_SIZE);
@ -601,11 +603,11 @@ double PairGranular::init_one(int i, int j)
double cutoff = 0.0;
if (setflag[i][j] == 0) {
if ((models[i][i]->normal_model.name != models[j][j]->normal_model.name) ||
(models[i][i]->damping_model.name != models[j][j]->damping_model.name) ||
(models[i][i]->tangential_model.name != models[j][j]->tangential_model.name) ||
(models[i][i]->rolling_model.name != models[j][j]->rolling_model.name) ||
(models[i][i]->twisting_model.name != models[j][j]->twisting_model.name)) {
if ((models[i][i]->normal_model->name != models[j][j]->normal_model->name) ||
(models[i][i]->damping_model->name != models[j][j]->damping_model->name) ||
(models[i][i]->tangential_model->name != models[j][j]->tangential_model->name) ||
(models[i][i]->rolling_model->name != models[j][j]->rolling_model->name) ||
(models[i][i]->twisting_model->name != models[j][j]->twisting_model->name)) {
error->all(FLERR,"Granular pair style functional forms are different, "
"cannot mix coefficients for types {} and {}. \n"
"This combination must be set explicitly via a "
@ -614,12 +616,12 @@ double PairGranular::init_one(int i, int j)
vec_models.push_back(ContactModel());
models[i][j] = models[j][i] = & vec_models.back();
vec_models.back().init(models[i][i]->normal_model.name, Contact::NORMAL);
vec_models.back().init(models[i][i]->tangential_model.name, Contact::TANGENTIAL);
vec_models.back().init(models[i][i]->damping_model.name, Contact::DAMPING);
vec_models.back().init(models[i][i]->rolling_model.name, Contact::ROLLING);
vec_models.back().init(models[i][i]->twisting_model.name, Contact::TWISTING);
vec_models.back().init(models[i][i]->heat_model.name, Contact::HEAT);
vec_models.back().init_model(models[i][i]->normal_model->name, NORMAL);
vec_models.back().init_model(models[i][i]->tangential_model->name, TANGENTIAL);
vec_models.back().init_model(models[i][i]->damping_model->name, DAMPING);
vec_models.back().init_model(models[i][i]->rolling_model->name, ROLLING);
vec_models.back().init_model(models[i][i]->twisting_model->name, TWISTING);
vec_models.back().init_model(models[i][i]->heat_model->name, HEAT);
vec_models.back().mix_coeffs(models[i][i], models[j][j]);
}
@ -675,10 +677,16 @@ double PairGranular::init_one(int i, int j)
// Copy global options
models[i][j]->dt = models[j][i]->dt = dt;
models[i][j]->normal_model.history_index = models[j][i]->normal_model.history_index = normal_history_index;
models[i][j]->tangential_model.history_index = models[j][i]->tangential_model.history_index = tangential_history_index;
models[i][j]->rolling_model.history_index = models[j][i]->rolling_model.history_index = rolling_history_index;
models[i][j]->twisting_model.history_index = models[j][i]->twisting_model.history_index = twisting_history_index;
models[i][j]->normal_model->history_index = models[j][i]->normal_model->history_index = normal_history_index;
models[i][j]->tangential_model->history_index = models[j][i]->tangential_model->history_index = tangential_history_index;
models[i][j]->rolling_model->history_index = models[j][i]->rolling_model->history_index = roll_history_index;
models[i][j]->twisting_model->history_index = models[j][i]->twisting_model->history_index = twist_history_index;
models[i][j]->size_history = models[j][i]->size_history = size_history;
models[i][j]->init(); // Calculates cumulative properties of sub models
models[j][i]->init();
if (models[i][j]->nondefault_history_transfer) nondefault_history_transfer = 1;
return cutoff;
}
@ -718,7 +726,7 @@ void PairGranular::read_restart(FILE *fp)
if (setflag[i][j]) {
vec_models.push_back(ContactModel());
models[i][j] = & vec_models.back();
models[i][j]->read_restart();
models[i][j]->read_restart(fp);
}
}
}
@ -729,6 +737,12 @@ void PairGranular::read_restart(FILE *fp)
void PairGranular::reset_dt()
{
dt = update->dt;
for (int i = 1; i <= atom->ntypes; i++) {
for (int j = i; j <= atom->ntypes; j++) {
models[i][j]->dt = dt;
}
}
}
/* ---------------------------------------------------------------------- */
@ -819,7 +833,7 @@ double PairGranular::single(int i, int j, int itype, int jtype,
history = &allhistory[size_history*neighprev];
}
double forces[3], torquesi[3], torquesj[3];
double *forces, *torquesi, *torquesj;
models[itype][jtype]->calculate_forces();
forces = models[itype][jtype]->forces;
torquesi = models[itype][jtype]->torquesi;
@ -893,11 +907,16 @@ double PairGranular::memory_usage()
only needed if any history entries i-j are not just negative of j-i entries
------------------------------------------------------------------------- */
void PairGranular::transfer_history(double* source, double* target, int source_type, int target_type)
void PairGranular::transfer_history(double* source, double* target, int itype, int jtype)
{
models[itype][jtype]->transfer_history();
for (int i = 0; i < size_history; i++){
target[i] = history_transfer_factors[i]*source[i];
if (models[itype][jtype]->nondefault_history_transfer) {
for (int i = 0; i < size_history; i++) {
target[i] = models[itype][jtype]->transfer_history_factor[i] * source [i];
}
} else {
for (int i = 0; i < size_history; i++) {
target[i] = -source[i];
}
}
}

6
src/GRANULAR/pair_granular.h
View File

@ -22,7 +22,7 @@ PairStyle(granular,PairGranular);
#include "contact.h"
#include "pair.h"
#include "vector.h"
#include <vector>
namespace LAMMPS_NS {
@ -65,11 +65,10 @@ class PairGranular : public Pair {
int nmax; // allocated size of mass_rigid
void allocate();
void transfer_history(double *, double *) override;
void transfer_history(double *, double *, int, int) override;
private:
int size_history;
int *history_transfer_factors;
int heat_flag;
// contact models
@ -81,6 +80,7 @@ class PairGranular : public Pair {
// indices of history entries
int normal_history_index;
int damping_history_index;
int tangential_history_index;
int roll_history_index;
int twist_history_index;

129
src/contact.h
View File

@ -14,78 +14,89 @@
#ifndef LMP_CONTACT_H
#define LMP_CONTACT_H
#include "pointers.h" // IWYU pragma: export
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};
enum ModelType {
NORMAL = 0,
DAMPING = 1,
TANGENTIAL = 2,
ROLLING = 3,
TWISTING = 4,
HEAT = 5
};
#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
#define EPSILON 1e-10
// forward declaration
class NormalModel;
class DampingModel;
class TangentialModel;
class RollingModel;
class TwistingModel;
class HeatModel;
class SubModel;
class ContactModel {
public:
ContactModel();
void reset_contact();
bool check_contact();
void prep_contact();
void calculate_forces(double *, double *, double *, double *);
double calculate_heat();
double pulloff_distance(double, double);
class ContactModel : protected Pointers {
public:
ContactModel();
~ContactModel();
void init();
bool check_contact();
void reset_contact();
void prep_contact();
void calculate_forces();
double calculate_heat();
double pulloff_distance(double, double);
int normal_model, damping_model, tangential_model;
int roll_model, twist_model;
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;
void init_model(std::string, ModelType);
double radi, radj, meff, dt, Ti, Tj;
double *xi, *xj, *vi, *vj, *omegai, *omegaj;
int history_update, roll_history_index, twist_history_index;
void mix_coeffs(ContactModel*, ContactModel*);
double fs[3], fr[3], ft[3], magtortwist;
void write_restart(FILE *);
void read_restart(FILE *);
private:
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;
// Sub models
NormalModel *normal_model;
DampingModel *damping_model;
TangentialModel *tangential_model;
RollingModel *rolling_model;
TwistingModel *twisting_model;
HeatModel *heat_model;
SubModel *sub_models[6]; // Need to resize if we add more model flavors
int prep_flag, check_flag;
int mindlin_rescale, mindlin_force;
// Extra options
int beyond_contact, limit_damping, history_update;
double cutoff_type;
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 *);
// History variables
int size_history, nondefault_history_transfer;
double *transfer_history_factor;
double *history;
};
// Contact properties/output
double *forces, *torquesi, *torquesj;
double radi, radj, meff, dt, Ti, Tj, area;
double Fntot, magtortwist;
double *xi, *xj, *vi, *vj, *omegai, *omegaj;
double fs[3], fr[3], ft[3];
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;
protected:
int prep_flag, check_flag;
int nmodels;
};
} // namespace Contact
} // namespace LAMMPS_NS
#endif

6
src/fix_neigh_history.cpp
View File

@ -337,6 +337,7 @@ void FixNeighHistory::pre_exchange_newton()
int *ilist,*jlist,*numneigh,**firstneigh;
int *allflags;
double *allvalues,*onevalues,*jvalues;
int *type = atom->type;
// NOTE: all operations until very end are with
// nlocal_neigh <= current nlocal and nall_neigh
@ -430,7 +431,7 @@ void FixNeighHistory::pre_exchange_newton()
partner[j][m] = tag[i];
jvalues = &valuepartner[j][dnum*m];
if (pair->nondefault_history_transfer)
pair->transfer_history(onevalues,jvalues);
pair->transfer_history(onevalues,jvalues,type[i],type[j]);
else for (n = 0; n < dnum; n++) jvalues[n] = -onevalues[n];
}
}
@ -469,6 +470,7 @@ void FixNeighHistory::pre_exchange_no_newton()
int *ilist,*jlist,*numneigh,**firstneigh;
int *allflags;
double *allvalues,*onevalues,*jvalues;
int *type = atom->type;
// NOTE: all operations until very end are with nlocal_neigh <= current nlocal
// because previous neigh list was built with nlocal_neigh
@ -544,7 +546,7 @@ void FixNeighHistory::pre_exchange_no_newton()
partner[j][m] = tag[i];
jvalues = &valuepartner[j][dnum*m];
if (pair->nondefault_history_transfer)
pair->transfer_history(onevalues, jvalues);
pair->transfer_history(onevalues, jvalues,type[i],type[j]);
else for (n = 0; n < dnum; n++) jvalues[n] = -onevalues[n];
}
}

2
src/pair.h
View File

@ -211,7 +211,7 @@ class Pair : protected Pointers {
virtual void min_xf_pointers(int, double **, double **) {}
virtual void min_xf_get(int) {}
virtual void min_x_set(int) {}
virtual void transfer_history(double *, double *) {}
virtual void transfer_history(double *, double *, int, int) {}
virtual double atom2cut(int) { return 0.0; }
virtual double radii2cut(double, double) { return 0.0; }