diff --git a/doc/src/pair_granular.rst b/doc/src/pair_granular.rst index 9804c90cef..6e17d604fd 100644 --- a/doc/src/pair_granular.rst +++ b/doc/src/pair_granular.rst @@ -89,7 +89,7 @@ and their required arguments are: 4. *dmt* : E, :math:`\eta_{n0}` (or :math:`e`), :math:`\nu`, :math:`\gamma` 5. *jkr* : E, :math:`\eta_{n0}` (or :math:`e`), :math:`\nu`, :math:`\gamma` 6. *mdr* : :math:`E`, :math:`\nu`, :math:`Y`, :math:`\Delta\gamma`, -:math:`\psi_b`, :math:`\eta_{n0}` + :math:`\psi_b`, :math:`\eta_{n0}` Here, :math:`k_n` is spring stiffness (with units that depend on model choice, see below); :math:`\eta_{n0}` is a damping prefactor (or, in its @@ -178,8 +178,9 @@ two-part series :ref:`Zunker and Kamrin Part I ` and :ref:`Zunker and Kamrin Part II `. Further development and demonstrations of its application to industrially relevant powder compaction processes are presented in :ref:`Zunker et al. `. -If you use the *mdr* normal model, is recommended you use the *mdr* damping model -described below. +If you use the *mdr* normal model, the default, and recommended damping model is +the *mdr* damping model described below. The other currently supported damping model +for the *mdr* normal model is *velocity*, all others will result in an error. The model requires the following inputs: @@ -203,8 +204,8 @@ The model requires the following inputs: triggered. Lower values of :math:`\psi_b` delay the onset of the bulk elastic response. - 6. *Coefficient of restitution* :math:`0 \le e \le 1` : The coefficient of - restitution is a tunable parameter that controls damping in the normal direction. + 6. *Damping coefficent* :math:`\eta_{n0} \ge 0` : The damping coefficient + is a tunable parameter that controls damping in the normal direction. .. note:: @@ -216,7 +217,7 @@ The *mdr* model produces a nonlinear force-displacement response, therefore the critical timestep :math:`\Delta t` depends on the inputs and level of deformation. As a conservative starting point the timestep can be assumed to be dictated by the bulk elastic response such that -:math:`\Delta t = 0.35\sqrt{m/k_\textrm{bulk}}`, where :math:`m` is the mass of +:math:`\Delta t = 0.08\sqrt{m/k_\textrm{bulk}}`, where :math:`m` is the mass of the smallest particle and :math:`k_\textrm{bulk} = \kappa R_\textrm{min}` is an effective stiffness related to the bulk elastic response. Here, :math:`\kappa = E/(3(1-2\nu))` is the bulk modulus and @@ -334,7 +335,8 @@ for the damping model currently supported are: *mdr* contact model is defined. If the *damping* keyword is not specified, the *viscoelastic* model is -used by default. +used by default. The exception is when the normal model is set +to *mdr* then the *mdr* damping model will be used by default. For *damping velocity*, the normal damping is simply equal to the user-specified damping coefficient in the *normal* model: @@ -344,7 +346,7 @@ user-specified damping coefficient in the *normal* model: \eta_n = \eta_{n0} Here, :math:`\eta_{n0}` is the damping coefficient specified for the normal -contact model, in units of *mass*\ /\ *time*\ . +contact model, in units of *mass*\ /\ *time*\ . For *damping mass_velocity*, the normal damping is given by: @@ -1082,8 +1084,8 @@ a bulk elastic response. Journal of the Mechanics and Physics of Solids, **(Zunker et al, 2025)** Zunker, W., Dunatunga, S., Thakur, S., Tang, P., & Kamrin, K. (2025). Experimentally validated DEM for large -deformation powder compaction: mechanically-derived contact model and -screening of non-physical contacts. +deformation powder compaction: Mechanically-derived contact model and +screening of non-physical contacts. Powder Technology, 120972. .. _Luding2008: diff --git a/examples/granular/in.tableting.200 b/examples/granular/in.tableting.200 index b9443cc36d..d1ecc28aa1 100644 --- a/examples/granular/in.tableting.200 +++ b/examples/granular/in.tableting.200 @@ -34,7 +34,7 @@ variable YieldStress equal 1.9e5 variable PoissonsRatio equal 0.4 variable SurfaceEnergy equal 2 variable SurfaceEnergyWall equal 0.0 -variable damp equal 1.0 +variable damp equal 0.2 variable psi_b equal 0.5 # linear_history = k_t, x_gammat, mu_s diff --git a/src/GRANULAR/gran_sub_mod_damping.cpp b/src/GRANULAR/gran_sub_mod_damping.cpp index a171f70654..13041b1a0b 100644 --- a/src/GRANULAR/gran_sub_mod_damping.cpp +++ b/src/GRANULAR/gran_sub_mod_damping.cpp @@ -77,7 +77,17 @@ GranSubModDampingVelocity::GranSubModDampingVelocity(GranularModel *gm, LAMMPS * double GranSubModDampingVelocity::calculate_forces() { - damp_prefactor = damp; + 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; + } return -damp_prefactor * gm->vnnr; } @@ -181,7 +191,6 @@ void GranSubModDampingCoeffRestitution::init() GranSubModDampingMDR::GranSubModDampingMDR(GranularModel *gm, LAMMPS *lmp) : GranSubModDamping(gm, lmp) { - contact_radius_flag = 1; } /* ---------------------------------------------------------------------- */ diff --git a/src/GRANULAR/gran_sub_mod_normal.cpp b/src/GRANULAR/gran_sub_mod_normal.cpp index 9eb30e2a0e..ec75fd2f71 100644 --- a/src/GRANULAR/gran_sub_mod_normal.cpp +++ b/src/GRANULAR/gran_sub_mod_normal.cpp @@ -73,8 +73,9 @@ static const char cite_mdr[] = " author = {Zunker, William and Dunatunga, Sachith and Thakur, Subhash and Tang, Pingjun and Kamrin, Ken},\n" " title = {Experimentally validated DEM for large deformation powder compaction:\n" " mechanically-derived contact model and screening of non-physical contacts},\n" + " journal = {Powder Technology},\n" " year = {2025},\n" - " journal = {engrXiv},\n" + " pages = {120972},\n" "}\n\n"; /* ---------------------------------------------------------------------- @@ -958,19 +959,18 @@ double GranSubModNormalMDR::calculate_forces() double damp_scale; if (gm->contact_type != PAIR) { a_damp = a_damp/2.0; - damp_scale = 0.5*sqrt(gm->meff * 2.0 * Eeff2particle * a_damp); + damp_scale = sqrt(gm->meff * 2.0 * Eeff2particle * a_damp); double *deltao_offset = &history[DELTAO_0]; const double wfm = std::exp(10.7*(*deltao_offset)/Rinitial[gm->i] - 10.0) + 1.0; // wall force magnifier - //const double wfm = 1.0; F = wij * F0 * wfm; } else { - damp_scale = 0.5*sqrt(gm->meff * 2.0 * Eeff * a_damp); + damp_scale = sqrt(gm->meff * 2.0 * Eeff * a_damp); F = wij * (F0 + F1) * 0.5; } if (history_update) { double *damp_scale_offset = & history[DAMP_SCALE]; - (a_damp < 0.0) ? *damp_scale_offset = 0.0 : *damp_scale_offset = damp_scale; + (a_damp <= 0.0) ? *damp_scale_offset = 0.0 : *damp_scale_offset = damp_scale; } return F;