diff --git a/doc/src/pair_granular.rst b/doc/src/pair_granular.rst index bc469412d9e..597869a05b4 100644 --- a/doc/src/pair_granular.rst +++ b/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) `, 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: diff --git a/doc/src/replicate.rst b/doc/src/replicate.rst index eb62fbfa210..8713ea2f95e 100644 --- a/doc/src/replicate.rst +++ b/doc/src/replicate.rst @@ -59,7 +59,7 @@ atoms. The bond discussion which follows only refers to models with permanent covalent bonds typically defined in LAMMPS via a data - file. It is not relevant to sytems modeled with many-body + file. It is not relevant to systems modeled with many-body potentials which can define bonds on-the-fly, based on the current positions of nearby atoms, e.g. models using the :doc:`AIREBO ` or :doc:`ReaxFF ` potentials. @@ -141,7 +141,7 @@ To work around this restriction two options are possible. (1) Fixes which use the stored data in the restart file can be defined before replication and then deleted via the :doc:`unfix ` command and re-defined after it. Or (2) the restart file can be converted to a -data file (which deletes the stored fix infomation) and fixes defined +data file (which deletes the stored fix information) and fixes defined after the replicate command. In both these scenarios, the per-atom fix information in the restart file is lost. diff --git a/doc/src/write_data.rst b/doc/src/write_data.rst index 516685f4fea..61fada70df9 100644 --- a/doc/src/write_data.rst +++ b/doc/src/write_data.rst @@ -12,14 +12,14 @@ Syntax * file = name of data file to write out * zero or more keyword/value pairs may be appended -* keyword = *pair* or *nocoeff* or *nofix* or *nolabelmap* +* keyword = *nocoeff* or *nofix* or *nolabelmap* or *triclinic/general* or *types* or *pair* .. parsed-literal:: *nocoeff* = do not write out force field info *nofix* = do not write out extra sections read by fixes *nolabelmap* = do not write out type labels - *triclinic/general = write data file in general triclinic format + *triclinic/general* = write data file in general triclinic format *types* value = *numeric* or *labels* *pair* value = *ii* or *ij* *ii* = write one line of pair coefficient info per atom type diff --git a/examples/granular/data.particles b/examples/granular/data.particles new file mode 100644 index 00000000000..c9f3bd7a9cf --- /dev/null +++ b/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 diff --git a/examples/granular/in.restitution b/examples/granular/in.restitution new file mode 100644 index 00000000000..e441ed67a75 --- /dev/null +++ b/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 diff --git a/examples/granular/log.13May23.restitution.g++.1 b/examples/granular/log.13May23.restitution.g++.1 new file mode 100644 index 00000000000..e51709d10dc --- /dev/null +++ b/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 diff --git a/src/EXTRA-FIX/fix_deform_pressure.cpp b/src/EXTRA-FIX/fix_deform_pressure.cpp index e2bcdd7f4e7..672f097c2da 100644 --- a/src/EXTRA-FIX/fix_deform_pressure.cpp +++ b/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) diff --git a/src/GRANULAR/gran_sub_mod_damping.cpp b/src/GRANULAR/gran_sub_mod_damping.cpp index 4386ed71fc7..df7146a6198 100644 --- a/src/GRANULAR/gran_sub_mod_damping.cpp +++ b/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 #include 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; +} diff --git a/src/GRANULAR/gran_sub_mod_damping.h b/src/GRANULAR/gran_sub_mod_damping.h index db5bb43ca5e..c931e385cca 100644 --- a/src/GRANULAR/gran_sub_mod_damping.h +++ b/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 diff --git a/src/GRANULAR/gran_sub_mod_normal.h b/src/GRANULAR/gran_sub_mod_normal.h index 6f2f3aabbe4..c1ed36b6e1c 100644 --- a/src/GRANULAR/gran_sub_mod_normal.h +++ b/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; diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index d5179a19b7d..cb19a8c0250 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/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]; diff --git a/src/ML-IAP/mliap_unified.cpp b/src/ML-IAP/mliap_unified.cpp index 7697204e441..3f0d0b2a8f0 100644 --- a/src/ML-IAP/mliap_unified.cpp +++ b/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; diff --git a/src/create_atoms.cpp b/src/create_atoms.cpp index 2912701159c..ded040b567d 100644 --- a/src/create_atoms.cpp +++ b/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);