From 210cddd94e2fefefa48dbdee3814ddc4ee7a79a5 Mon Sep 17 00:00:00 2001 From: William Zunker Date: Tue, 8 Apr 2025 18:32:45 -0400 Subject: [PATCH] working on adding damping types to mdr damping --- src/GRANULAR/gran_sub_mod_damping.cpp | 27 +++++++++++++++------------ src/GRANULAR/gran_sub_mod_damping.h | 1 + src/GRANULAR/gran_sub_mod_normal.cpp | 16 +++++++++------- src/GRANULAR/gran_sub_mod_normal.h | 2 ++ 4 files changed, 27 insertions(+), 19 deletions(-) diff --git a/src/GRANULAR/gran_sub_mod_damping.cpp b/src/GRANULAR/gran_sub_mod_damping.cpp index da9aa81079..a5dc5745a4 100644 --- a/src/GRANULAR/gran_sub_mod_damping.cpp +++ b/src/GRANULAR/gran_sub_mod_damping.cpp @@ -45,6 +45,11 @@ GranSubModDamping::GranSubModDamping(GranularModel *gm, LAMMPS *lmp) : GranSubMo void GranSubModDamping::init() { damp = gm->normal_model->get_damp(); + if (gm->normal_model->name == "mdr") { + if (gm->tangential_model->name != "mdr") + error->all(FLERR, "Only damping mdr may be used with the mdr normal model"); + damp_type = gm->normal_model->get_damp_type(); + } } /* ---------------------------------------------------------------------- @@ -77,17 +82,7 @@ GranSubModDampingVelocity::GranSubModDampingVelocity(GranularModel *gm, LAMMPS * double GranSubModDampingVelocity::calculate_forces() { - if (gm->normal_model->name == "mdr") { - using namespace Granular_MDR_NS; - double *history = & gm->history[gm->normal_model->history_index]; - if (history[DAMP_SCALE] == 0.0) { - damp_prefactor == 0.0; - } else { - damp_prefactor = damp; - } - } else { - damp_prefactor = damp; - } + damp_prefactor = damp; return -damp_prefactor * gm->vnnr; } @@ -209,6 +204,14 @@ double GranSubModDampingMDR::calculate_forces() { using namespace Granular_MDR_NS; double *history = & gm->history[gm->normal_model->history_index]; - damp_prefactor = damp * history[DAMP_SCALE]; + if (damp_type == 1) { + damp_prefactor = damp * history[DAMP_SCALE]; + } else if (damp_type == 2) { + if (history[DAMP_SCALE] == 0.0) { + damp_prefactor == 0.0; + } else { + damp_prefactor = damp; + } + } return -damp_prefactor * gm->vnnr; } diff --git a/src/GRANULAR/gran_sub_mod_damping.h b/src/GRANULAR/gran_sub_mod_damping.h index cd7ee3bf1a..c73a4a8e5e 100644 --- a/src/GRANULAR/gran_sub_mod_damping.h +++ b/src/GRANULAR/gran_sub_mod_damping.h @@ -42,6 +42,7 @@ namespace Granular_NS { protected: double damp_prefactor; double damp; + int damp_type; }; /* ---------------------------------------------------------------------- */ diff --git a/src/GRANULAR/gran_sub_mod_normal.cpp b/src/GRANULAR/gran_sub_mod_normal.cpp index ec75fd2f71..ecf53ea24c 100644 --- a/src/GRANULAR/gran_sub_mod_normal.cpp +++ b/src/GRANULAR/gran_sub_mod_normal.cpp @@ -441,7 +441,7 @@ GranSubModNormalMDR::GranSubModNormalMDR(GranularModel *gm, LAMMPS *lmp) : { if (lmp->citeme) lmp->citeme->add(cite_mdr); - num_coeffs = 6; + num_coeffs = 7; contact_radius_flag = 1; size_history = 27; nsvector = 1; @@ -468,12 +468,13 @@ GranSubModNormalMDR::~GranSubModNormalMDR() void GranSubModNormalMDR::coeffs_to_local() { - E = coeffs[0]; // Young's modulus - nu = coeffs[1]; // Poisson's ratio - Y = coeffs[2]; // yield stress - gamma = coeffs[3]; // effective surface energy - psi_b = coeffs[4]; // bulk response trigger based on ratio of remaining free area: A_{free}/A_{total} - damp = coeffs[5]; // coefficent of restitution + E = coeffs[0]; // Young's modulus + nu = coeffs[1]; // Poisson's ratio + Y = coeffs[2]; // yield stress + gamma = coeffs[3]; // effective surface energy + psi_b = coeffs[4]; // bulk response trigger based on ratio of remaining free area: A_{free}/A_{total} + damp = coeffs[5]; // coefficent of restitution + damp_type = coeffs[6]; // damping type 1 = mdr or 2 = velocity if (E <= 0.0) error->all(FLERR, "Illegal MDR normal model, Young's modulus must be greater than 0"); if (nu < 0.0 || nu > 0.5) error->all(FLERR, "Illegal MDR normal model, Poisson's ratio must be between 0 and 0.5"); @@ -481,6 +482,7 @@ void GranSubModNormalMDR::coeffs_to_local() if (gamma < 0.0) error->all(FLERR, "Illegal MDR normal model, effective surface energy must be greater than or equal to 0"); if (psi_b < 0.0 || psi_b > 1.0) error->all(FLERR, "Illegal MDR normal model, psi_b must be between 0 and 1.0"); if (damp < 0.0) error->all(FLERR, "Illegal MDR normal model, damping coefficent must be greater than or equal to 0"); + if (damp_type != 1 && damp_type != 2) error->all(FLERR, "Illegal MDR normal model, damping type must an integer equal to 1 (mdr) or 2 (velocity)"); G = E / (2.0 * (1.0 + nu)); // shear modulus kappa = E / (3.0 * (1.0 - 2.0 * nu)); // bulk modulus diff --git a/src/GRANULAR/gran_sub_mod_normal.h b/src/GRANULAR/gran_sub_mod_normal.h index 8142d88cea..0b7bb664d5 100644 --- a/src/GRANULAR/gran_sub_mod_normal.h +++ b/src/GRANULAR/gran_sub_mod_normal.h @@ -41,6 +41,7 @@ namespace Granular_NS { int get_cohesive_flag() const { return cohesive_flag; } double get_damp() const { return damp; } + double get_damp_type() const { return damp_type; } double get_emod() const { return Emod; } double get_fncrit() const { return Fncrit; } int get_material_properties() const { return material_properties; } @@ -51,6 +52,7 @@ namespace Granular_NS { protected: double damp; // argument historically needed by damping // typically (but not always) equals eta_n0 + int damp_type; // damping type is only used by normal mdr model double Emod, poiss; double Fncrit; int material_properties, cohesive_flag;