Unifying desc/code for tsuji and coeff restitution

doc/src/pair_granular.rst

@@ -253,26 +253,26 @@ of the normal contact model parameters should be between 0 and 1, but
 no error check is performed on this.
 
 The *coeff_restitution* model is useful when a specific normal coefficient of
-restitution :math:`e` is required. In these models, the normal coefficient of
+restitution :math:`e` is required. It operates much like the *Tsuji* model
-restitution :math:`e` is specified as an input in place of the usual :math:`\eta_{n0}`
+but, the normal coefficient of restitution :math:`e` is specified as an input
-value in the normal model. Following the approach of
+in place of the usual :math:`\eta_{n0}` value in the normal model. Following
-:ref:`(Brilliantov et al) <Brill1996>`, when using the *hooke* normal model,
+the approach of :ref:`(Brilliantov et al) <Brill1996>`, when using the *hooke*
-*coeff_restitution* then calculates the damping coefficient as:
+normal model, *coeff_restitution* then calculates the damping coefficient as:
 
 .. math::
 
-   \eta_n = \sqrt{\frac{4m_{eff}k_n}{1+\left( \frac{\pi}{\log(e)}\right)^2}} ,
+   \eta_n = \sqrt{\frac{4m_{eff}k_{nd}}{1+\left( \frac{\pi}{\log(e)}\right)^2}} ,
 
+where :math:`k_{nd}` is the same stiffness defined in the above *Tsuji* model.
 For any other normal model, e.g. the *hertz* and *hertz/material* models, the damping
 coefficient is:
 
 .. math::
 
-   \eta_n = -2\sqrt{\frac{5}{6}}\frac{\log(e)}{\sqrt{\pi^2+(\log(e))^2}}(R_{eff} \delta_{ij})^{\frac{1}{4}}\sqrt{\frac{3}{2}k_n m_{eff}} ,
+   \eta_n = -2\sqrt{\frac{5}{6}}\frac{\log(e)}{\sqrt{\pi^2+(\log(e))^2}}(R_{eff} \delta_{ij})^{\frac{1}{4}}\sqrt{\frac{3}{2}k_{nd} m_{eff}} ,
 
-where :math:`k_n = \frac{4}{3} E_{eff}` for the *hertz/material* model. Since
+Since *coeff_restitution* accounts for the effective mass, effective radius,
-*coeff_restitution* accounts for the effective mass, effective radius, and
+and pairwise overlaps (except when used with the *hooke* normal model) when calculating
-pairwise overlaps (except when used with the *hooke* normal model) when calculating
 the damping coefficient, it accurately reproduces the specified coefficient of
 restitution for both monodisperse and polydisperse particle pairs.  This damping
 model is not compatible with cohesive normal models such as *JKR* or *DMT*.

src/GRANULAR/gran_sub_mod_damping.cpp

@@ -152,9 +152,8 @@ double GranSubModDampingTsuji::calculate_forces()
 ------------------------------------------------------------------------- */
 
 GranSubModDampingCoeffRestitution::GranSubModDampingCoeffRestitution(GranularModel *gm, LAMMPS *lmp) :
-    GranSubModDamping(gm, lmp)
+    GranSubModDampingTsuji(gm, lmp)
 {
-  allow_cohesion = 0;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -171,17 +170,3 @@ void GranSubModDampingCoeffRestitution::init()
     damp /= sqrt(MY_PI * MY_PI + logcor * logcor);
   }
 }
-
-/* ---------------------------------------------------------------------- */
-
-double GranSubModDampingCoeffRestitution::calculate_forces()
-{
-  // in case argument <= 0 due to precision issues
-  double sqrt1;
-  if (gm->delta > 0.0)
-    sqrt1 = MAX(0.0, gm->meff * gm->Fnormal / gm->delta);
-  else
-    sqrt1 = 0.0;
-  damp_prefactor = damp * sqrt(sqrt1);
-  return -damp_prefactor * gm->vnnr;
-}

src/GRANULAR/gran_sub_mod_damping.h

@@ -87,11 +87,10 @@ namespace Granular_NS {
 
   /* ---------------------------------------------------------------------- */
 
-  class GranSubModDampingCoeffRestitution : public GranSubModDamping {
+  class GranSubModDampingCoeffRestitution : public GranSubModDampingTsuji {
    public:
     GranSubModDampingCoeffRestitution(class GranularModel *, class LAMMPS *);
     void init() override;
-    double calculate_forces() override;
   };
 
   /* ---------------------------------------------------------------------- */