cleaned mdr damping and updated velocity damping for mdr use

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 <Zunker2024I>` and
 :ref:`Zunker and Kamrin Part II <Zunker2024II>`. Further development
 and demonstrations of its application to industrially relevant powder
 compaction processes are presented in :ref:`Zunker et al. <Zunker2025>`.
-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:
 

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

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;
 }
 
 /* ---------------------------------------------------------------------- */

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;