Merge pull request #4068 from dhairyaiitb/develop

Coefficient of restitution based damping in granular models

doc/src/pair_granular.rst

@@ -187,6 +187,7 @@ for the damping model currently supported are:
 2. *mass_velocity*
 3. *viscoelastic*
 4. *tsuji*
+5. *coeff_restitution*
 
 If the *damping* keyword is not specified, the *viscoelastic* model is
 used by default.
@@ -248,6 +249,29 @@ The dimensionless coefficient of restitution :math:`e` specified as part
 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 specified as an input. Following the approach of
+:ref:`(Brilliantov et al) <Brill1996>`, when using the *hooke* normal model,
+*coeff_restitution* calculates the damping coefficient as:
+
+.. math::
+
+   \eta_n = \sqrt{\frac{4m_{eff}k_n}{1+\left( \frac{\pi}{\log(e)}\right)^2}} ,
+
+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}} ,
+
+where :math:`k_n = \frac{4}{3} E_{eff}` for the *hertz/material* model. Since
+*coeff_restitution* accounts for the effective mass, effective radius, and
+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.
+
 The total normal force is computed as the sum of the elastic and
 damping components:
 

examples/granular/data.particles

@@ -0,0 +1,18 @@
+Python generated LAMMPS data file
+
+2 atoms
+1 atom types
+
+0 0.08 xlo xhi
+0 0.04 ylo yhi
+0 0.08 zlo zhi
+
+Atoms # sphere
+
+1 1 0.004 2500 0.04 0.02 0.04 0 0 0
+2 1 0.004 2500 0.04 0.02 0.04416 0 0 0
+
+Velocities
+
+1 0.0  0.0 1 0 0 0
+2 0.0  0.0 -1 0 0 0

examples/granular/in.restitution

@@ -0,0 +1,24 @@
+units	        si
+atom_style	    sphere
+
+boundary    	p p f
+region		    box block 0 80e-3 0 40e-3 0 80e-3 open 3 open 4
+create_box	    2 box
+
+read_data	    data.particles add append
+group 		    mb type 1
+
+pair_style granular
+pair_coeff * * hertz/material 1e6 0.8 0.4 tangential mindlin  NULL 0.0 0.5 damping coeff_restitution
+# pair_coeff * * hooke 1e6 0.5 tangential mindlin  1 1.0 0.0 damping coeff_restitution
+comm_modify     vel yes
+
+timestep 	    1e-9
+fix		        1 all nve/sphere
+compute		    s all stress/atom NULL pair
+
+#dump		    1 all custom 2000000 op.dump id x y z vx vy vz
+#dump_modify     1 pad 8
+thermo_style	custom step ke
+run_style 	    verlet
+run		        10000000

examples/granular/log.13May23.restitution.g++.1

@@ -0,0 +1,80 @@
+LAMMPS (17 Apr 2024 - Development - patch_17Apr2024-93-g4e7bddaa0b)
+  using 1 OpenMP thread(s) per MPI task
+units	        si
+atom_style	    sphere
+
+boundary    	p p f
+region		    box block 0 80e-3 0 40e-3 0 80e-3 open 3 open 4
+create_box	    2 box
+Created orthogonal box = (0 0 0) to (0.08 0.04 0.08)
+  1 by 1 by 1 MPI processor grid
+
+read_data	    data.particles add append
+Reading data file ...
+  orthogonal box = (0 0 0) to (0.08 0.04 0.08)
+  1 by 1 by 1 MPI processor grid
+  reading atoms ...
+  2 atoms
+  reading velocities ...
+  2 velocities
+  read_data CPU = 0.002 seconds
+group 		    mb type 1
+2 atoms in group mb
+
+pair_style granular
+pair_coeff * * hertz/material 1e6 0.8 0.4 tangential mindlin  NULL 0.0 0.5 damping coeff_restitution
+# pair_coeff * * hooke 1e6 0.5 tangential mindlin  1 1.0 0.0 damping coeff_restitution
+comm_modify     vel yes
+
+timestep 	    1e-9
+fix		        1 all nve/sphere
+compute		    s all stress/atom NULL pair
+
+#dump		    1 all custom 2000000 op.dump id x y z vx vy vz
+#dump_modify     1 pad 8
+thermo_style	custom step ke
+run_style 	    verlet
+run		        10000000
+Generated 0 of 1 mixed pair_coeff terms from geometric mixing rule
+Neighbor list info ...
+  update: every = 1 steps, delay = 0 steps, check = yes
+  max neighbors/atom: 2000, page size: 100000
+  master list distance cutoff = 0.005
+  ghost atom cutoff = 0.005
+  binsize = 0.0025, bins = 32 16 32
+  1 neighbor lists, perpetual/occasional/extra = 1 0 0
+  (1) pair granular, perpetual
+      attributes: half, newton on, size, history
+      pair build: half/size/bin/atomonly/newton
+      stencil: half/bin/3d
+      bin: standard
+Per MPI rank memory allocation (min/avg/max) = 10.1 | 10.1 | 10.1 Mbytes
+   Step         KinEng    
+         0   8.3775804e-05
+  10000000   5.3616513e-05
+Loop time of 5.99782 on 1 procs for 10000000 steps with 2 atoms
+
+77.9% CPU use with 1 MPI tasks x 1 OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 0.60235    | 0.60235    | 0.60235    |   0.0 | 10.04
+Neigh   | 0.00021965 | 0.00021965 | 0.00021965 |   0.0 |  0.00
+Comm    | 1.7939     | 1.7939     | 1.7939     |   0.0 | 29.91
+Output  | 2.5955e-05 | 2.5955e-05 | 2.5955e-05 |   0.0 |  0.00
+Modify  | 1.7622     | 1.7622     | 1.7622     |   0.0 | 29.38
+Other   |            | 1.839      |            |       | 30.66
+
+Nlocal:              2 ave           2 max           2 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Nghost:              0 ave           0 max           0 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Neighs:              0 ave           0 max           0 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+
+Total # of neighbors = 0
+Ave neighs/atom = 0
+Neighbor list builds = 14
+Dangerous builds = 0
+Total wall time: 0:00:06

src/EXTRA-FIX/fix_deform_pressure.cpp

@@ -438,8 +438,6 @@ void FixDeformPressure::init()
   // error if style PRESSURE/ERATEER for yz, can't calculate if box flip occurs
 
   if (set[3].style && set[5].style) {
-    int flag = 0;
-    double lo,hi;
     if (flipflag && set[3].style == PRESSURE)
       error->all(FLERR, "Fix {} cannot use yz pressure with xy", style);
     if (flipflag && set[3].style == ERATERS)

src/GRANULAR/gran_sub_mod_damping.cpp

@@ -16,16 +16,23 @@
 #include "gran_sub_mod_normal.h"
 #include "granular_model.h"
 #include "math_special.h"
+#include "math_const.h"
+
+#include <cmath>
 
 #include <cmath>
 
 using namespace LAMMPS_NS;
 using namespace Granular_NS;
+using namespace MathConst;
 
 using MathSpecial::cube;
 using MathSpecial::powint;
 using MathSpecial::square;
 
+static constexpr double TWOROOTFIVEBYSIX = 1.82574185835055380345;      // 2/sqrt(5/6)
+static constexpr double ROOTTHREEBYTWO = 1.22474487139158894067;        // sqrt(3/2)
+
 /* ----------------------------------------------------------------------
    Default damping model
 ------------------------------------------------------------------------- */
@@ -141,3 +148,31 @@ double GranSubModDampingTsuji::calculate_forces()
   damp_prefactor = damp * sqrt(sqrt1);
   return -damp_prefactor * gm->vnnr;
 }
+
+/* ----------------------------------------------------------------------
+   Coefficient of restitution damping
+------------------------------------------------------------------------- */
+
+GranSubModDampingCoeffRestitution::GranSubModDampingCoeffRestitution(GranularModel *gm, LAMMPS *lmp) :
+    GranSubModDamping(gm, lmp)
+{
+}
+
+void GranSubModDampingCoeffRestitution::init()
+{
+  // Calculate prefactor, assume Hertzian as default
+  double cor = gm->normal_model->get_damp();
+  double logcor = log(cor);
+  if (gm->normal_model->name == "hooke") {
+    damp = -2 * logcor / sqrt(MY_PI * MY_PI + logcor * logcor);
+  } else {
+    damp = -ROOTTHREEBYTWO * TWOROOTFIVEBYSIX * logcor;
+    damp /= sqrt(MY_PI * MY_PI + logcor * logcor);
+  }
+}
+
+double GranSubModDampingCoeffRestitution::calculate_forces()
+{
+  damp_prefactor = damp * sqrt(gm->meff * gm->Fnormal / gm->delta);
+  return -damp_prefactor * gm->vnnr;
+}

src/GRANULAR/gran_sub_mod_damping.h

@@ -18,6 +18,7 @@ GranSubModStyle(velocity,GranSubModDampingVelocity,DAMPING);
 GranSubModStyle(mass_velocity,GranSubModDampingMassVelocity,DAMPING);
 GranSubModStyle(viscoelastic,GranSubModDampingViscoelastic,DAMPING);
 GranSubModStyle(tsuji,GranSubModDampingTsuji,DAMPING);
+GranSubModStyle(coeff_restitution,GranSubModDampingCoeffRestitution,DAMPING);
 // clang-format on
 #else
 
@@ -84,6 +85,17 @@ namespace Granular_NS {
     double calculate_forces() override;
   };
 
+  /* ---------------------------------------------------------------------- */
+
+  class GranSubModDampingCoeffRestitution : public GranSubModDamping {
+   public:
+    GranSubModDampingCoeffRestitution(class GranularModel *, class LAMMPS *);
+    void init() override;
+    double calculate_forces() override;
+  };
+
+  /* ---------------------------------------------------------------------- */
+
 }    // namespace Granular_NS
 }    // namespace LAMMPS_NS
 

src/GRANULAR/gran_sub_mod_normal.h

@@ -38,7 +38,7 @@ namespace Granular_NS {
     virtual double calculate_contact_radius();
     virtual double calculate_forces() = 0;
 
-    double get_cohesive_flag() const { return cohesive_flag; }
+    int get_cohesive_flag() const { return cohesive_flag; }
     double get_damp() const { return damp; }
     double get_emod() const { return Emod; }
     double get_fncrit() const { return Fncrit; }
@@ -49,6 +49,7 @@ namespace Granular_NS {
 
    protected:
     double damp;    // argument historically needed by damping
+                    // typically (but not always) equals eta_n0
     double Emod, poiss;
     double Fncrit;
     int material_properties, cohesive_flag;

src/GRANULAR/pair_granular.cpp

@@ -769,7 +769,7 @@ double PairGranular::single(int i, int j, int itype, int jtype,
 
   // apply forces & torques
   // Calculate normal component, normalized by r
-  fforce = model->Fnormal * model->rinv;
+  fforce = model->Fntot * model->rinv;
 
   // set single_extra quantities
   svector[0] = model->fs[0];

src/ML-IAP/mliap_unified.cpp

@@ -246,7 +246,6 @@ void LAMMPS_NS::update_pair_energy(MLIAPData *data, double *eij)
 {
   double e_total = 0.0;
   const auto nlistatoms = data->nlistatoms;
-  const auto nlocal = data->nlocal;
   for (int ii = 0; ii < nlistatoms; ii++) data->eatoms[ii] = 0;
 
   for (int ii = 0; ii < data->npairs; ii++) {
@@ -268,7 +267,6 @@ void LAMMPS_NS::update_pair_forces(MLIAPData *data, double *fij)
 {
   //Bugfix: need to account for Null atoms in local atoms
   //const auto nlistatoms = data->nlistatoms;
-  const auto nlocal = data->nlocal;
   double **f = data->f;
   for (int ii = 0; ii < data->npairs; ii++) {
     int ii3 = ii * 3;

src/create_atoms.cpp

@@ -1505,7 +1505,7 @@ void CreateAtoms::get_xmol(double *center)
 
   onemol->quat_external = quatone;
 
-  int n, natoms = onemol->natoms;
+  int natoms = onemol->natoms;
   double xnew[3];
   for (int m = 0; m < natoms; m++) {
     MathExtra::matvec(rotmat, onemol->dx[m], xnew);