From 63f33aa3a428e81e2e5747cc3f9eb96182788c3b Mon Sep 17 00:00:00 2001 From: Dhairya Date: Sun, 4 Feb 2024 16:29:19 -0600 Subject: [PATCH 01/13] first --- src/GRANULAR/granular_model.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/GRANULAR/granular_model.cpp b/src/GRANULAR/granular_model.cpp index c1ad692fb3f..af095d1efcb 100644 --- a/src/GRANULAR/granular_model.cpp +++ b/src/GRANULAR/granular_model.cpp @@ -10,7 +10,7 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ - + /* ---------------------------------------------------------------------- This class contains a series of tools for DEM contacts Multiple models can be defined and used to calculate forces From 5137e8697250650f10076495bdf8d43e69d9b1b3 Mon Sep 17 00:00:00 2001 From: Dhairya Date: Mon, 5 Feb 2024 16:30:44 -0600 Subject: [PATCH 02/13] en models incorporated --- src/GRANULAR/gran_sub_mod_damping.cpp | 50 +++++++++++++++++++++++++++ src/GRANULAR/gran_sub_mod_damping.h | 22 ++++++++++++ 2 files changed, 72 insertions(+) diff --git a/src/GRANULAR/gran_sub_mod_damping.cpp b/src/GRANULAR/gran_sub_mod_damping.cpp index 7d6a02b8f0b..2ea9eed0aa2 100644 --- a/src/GRANULAR/gran_sub_mod_damping.cpp +++ b/src/GRANULAR/gran_sub_mod_damping.cpp @@ -17,6 +17,8 @@ #include "granular_model.h" #include "math_special.h" +#include + using namespace LAMMPS_NS; using namespace Granular_NS; @@ -24,6 +26,10 @@ using MathSpecial::cube; using MathSpecial::powint; using MathSpecial::square; +static constexpr double PISQ = 9.86960440108935799230; +static constexpr double TWOROOTFIVEBYSIX = 1.82574185835055380345; +static constexpr double ROOTTHREEBYTWO = 1.22474487139158894067; + /* ---------------------------------------------------------------------- Default damping model ------------------------------------------------------------------------- */ @@ -133,3 +139,47 @@ double GranSubModDampingTsuji::calculate_forces() damp_prefactor = damp * sqrt(gm->meff * gm->Fnormal / gm->delta); return -damp_prefactor * gm->vnnr; } + +/* ---------------------------------------------------------------------- + enhooke damping +------------------------------------------------------------------------- */ + +GranSubModDampingEnHooke::GranSubModDampingEnHooke(GranularModel *gm, LAMMPS *lmp) : + GranSubModDamping(gm, lmp) +{ +} + +void GranSubModDampingEnHooke::init() +{ + double cor = gm->normal_model->get_damp(); + double logcor = log(cor); + damp = -2*logcor/sqrt(PISQ + logcor*logcor); +} + +double GranSubModDampingEnHooke::calculate_forces() +{ + damp_prefactor = damp * sqrt(gm->meff * gm->Fnormal / gm->delta); + return -damp_prefactor * gm->vnnr; +} + +/* ---------------------------------------------------------------------- + enhertz damping +------------------------------------------------------------------------- */ + +GranSubModDampingEnHertz::GranSubModDampingEnHertz(GranularModel *gm, LAMMPS *lmp) : + GranSubModDamping(gm, lmp) +{ +} + +void GranSubModDampingEnHertz::init() +{ + double cor = gm->normal_model->get_damp(); + double logcor = log(cor); + damp = -ROOTTHREEBYTWO*TWOROOTFIVEBYSIX*logcor/sqrt(PISQ + logcor*logcor); +} + +double GranSubModDampingEnHertz::calculate_forces() +{ + damp_prefactor = damp * sqrt(gm->meff * gm->Fnormal / gm->delta); + return -damp_prefactor * gm->vnnr; +} \ No newline at end of file diff --git a/src/GRANULAR/gran_sub_mod_damping.h b/src/GRANULAR/gran_sub_mod_damping.h index db5bb43ca5e..33e88718fd3 100644 --- a/src/GRANULAR/gran_sub_mod_damping.h +++ b/src/GRANULAR/gran_sub_mod_damping.h @@ -18,6 +18,8 @@ GranSubModStyle(velocity,GranSubModDampingVelocity,DAMPING); GranSubModStyle(mass_velocity,GranSubModDampingMassVelocity,DAMPING); GranSubModStyle(viscoelastic,GranSubModDampingViscoelastic,DAMPING); GranSubModStyle(tsuji,GranSubModDampingTsuji,DAMPING); +GranSubModStyle(enhooke,GranSubModDampingEnHooke,DAMPING); +GranSubModStyle(enhertz,GranSubModDampingEnHertz,DAMPING); // clang-format on #else @@ -84,6 +86,26 @@ namespace Granular_NS { double calculate_forces() override; }; + /* ---------------------------------------------------------------------- */ + + class GranSubModDampingEnHooke : public GranSubModDamping { + public: + GranSubModDampingEnHooke(class GranularModel *, class LAMMPS *); + void init() override; + double calculate_forces() override; + }; + + /* ---------------------------------------------------------------------- */ + + class GranSubModDampingEnHertz : public GranSubModDamping { + public: + GranSubModDampingEnHertz(class GranularModel *, class LAMMPS *); + void init() override; + double calculate_forces() override; + }; + + /* ---------------------------------------------------------------------- */ + } // namespace Granular_NS } // namespace LAMMPS_NS From febb3671d82c499dbe9a381a0a65d30572c409b7 Mon Sep 17 00:00:00 2001 From: Dhairya Date: Tue, 6 Feb 2024 18:03:17 -0600 Subject: [PATCH 03/13] removed whitespace --- src/GRANULAR/granular_model.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/GRANULAR/granular_model.cpp b/src/GRANULAR/granular_model.cpp index af095d1efcb..c1ad692fb3f 100644 --- a/src/GRANULAR/granular_model.cpp +++ b/src/GRANULAR/granular_model.cpp @@ -10,7 +10,7 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ - + /* ---------------------------------------------------------------------- This class contains a series of tools for DEM contacts Multiple models can be defined and used to calculate forces From 26cff47386d94bdda1ff814ddb4f247120ce97be Mon Sep 17 00:00:00 2001 From: Dhairya Date: Wed, 7 Feb 2024 09:45:29 -0600 Subject: [PATCH 04/13] Removed whitespace --- src/GRANULAR/gran_sub_mod_damping.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/GRANULAR/gran_sub_mod_damping.cpp b/src/GRANULAR/gran_sub_mod_damping.cpp index 2ea9eed0aa2..e57972db767 100644 --- a/src/GRANULAR/gran_sub_mod_damping.cpp +++ b/src/GRANULAR/gran_sub_mod_damping.cpp @@ -26,7 +26,7 @@ using MathSpecial::cube; using MathSpecial::powint; using MathSpecial::square; -static constexpr double PISQ = 9.86960440108935799230; +static constexpr double PISQ = 9.86960440108935799230; static constexpr double TWOROOTFIVEBYSIX = 1.82574185835055380345; static constexpr double ROOTTHREEBYTWO = 1.22474487139158894067; @@ -182,4 +182,4 @@ double GranSubModDampingEnHertz::calculate_forces() { damp_prefactor = damp * sqrt(gm->meff * gm->Fnormal / gm->delta); return -damp_prefactor * gm->vnnr; -} \ No newline at end of file +} From 194b45b7295bb9ad17c29dae482cdc51f1c575f9 Mon Sep 17 00:00:00 2001 From: Dhairya Date: Wed, 7 Feb 2024 09:52:02 -0600 Subject: [PATCH 05/13] Example file --- examples/granular/en_example/particles.dat | 18 ++++++++++++++++ examples/granular/en_example/start.lammps | 25 ++++++++++++++++++++++ 2 files changed, 43 insertions(+) create mode 100644 examples/granular/en_example/particles.dat create mode 100644 examples/granular/en_example/start.lammps diff --git a/examples/granular/en_example/particles.dat b/examples/granular/en_example/particles.dat new file mode 100644 index 00000000000..c9f3bd7a9cf --- /dev/null +++ b/examples/granular/en_example/particles.dat @@ -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/en_example/start.lammps b/examples/granular/en_example/start.lammps new file mode 100644 index 00000000000..2bf0f54ddc0 --- /dev/null +++ b/examples/granular/en_example/start.lammps @@ -0,0 +1,25 @@ +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 particles.dat 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 enhertz +# pair_coeff * * hertz 1e6 0.3 tangential mindlin 1e4 1.0 0.5 damping enhertz +# pair_coeff * * hooke 1e6 0.5 tangential mindlin 1 1.0 0.0 damping enhooke +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 From f9ecdb5b546c1243b2669600cf2c2a25dec1b1fd Mon Sep 17 00:00:00 2001 From: Dhairya Date: Fri, 23 Feb 2024 15:21:56 -0600 Subject: [PATCH 06/13] Updated documentation --- doc/src/pair_granular.rst | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/doc/src/pair_granular.rst b/doc/src/pair_granular.rst index bc469412d9e..f0eb139cb72 100644 --- a/doc/src/pair_granular.rst +++ b/doc/src/pair_granular.rst @@ -187,6 +187,8 @@ for the damping model currently supported are: 2. *mass_velocity* 3. *viscoelastic* 4. *tsuji* +5. *enhooke* +6. *enhertz* If the *damping* keyword is not specified, the *viscoelastic* model is used by default. @@ -248,6 +250,21 @@ 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. +*enhooke* and *enhertz* models are useful for cases where 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) `, *enhooke* calculates the damping coefficient for the *hooke* model as: + +.. math:: + + \eta_n = \sqrt{\frac{4m_{eff}k_n}{1+\left( \frac{\pi}{\log(e)}\right)^2}}, + + +*enhertz* calculates the damping coefficient for the *hertz* and *hertz/material* models using: + +.. 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 these models calculate the damping coefficients by accounting for the effective mass, effective radius and pairwise overlaps (for *enhertz*), they accurately reproduce 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: From 84b6c6a088a2f36b80a7ebf6eb876c28120e1100 Mon Sep 17 00:00:00 2001 From: Dhairya Date: Fri, 5 Apr 2024 09:58:33 -0500 Subject: [PATCH 07/13] Added prefactors and errors for incorrect combinations. --- doc/src/pair_granular.rst | 10 +++++----- src/GRANULAR/gran_sub_mod_damping.cpp | 16 ++++++++-------- src/GRANULAR/gran_sub_mod_damping.h | 12 ++++++------ src/GRANULAR/granular_model.cpp | 4 ++++ 4 files changed, 23 insertions(+), 19 deletions(-) diff --git a/doc/src/pair_granular.rst b/doc/src/pair_granular.rst index f0eb139cb72..9eb34c9de84 100644 --- a/doc/src/pair_granular.rst +++ b/doc/src/pair_granular.rst @@ -187,8 +187,8 @@ for the damping model currently supported are: 2. *mass_velocity* 3. *viscoelastic* 4. *tsuji* -5. *enhooke* -6. *enhertz* +5. *hooke/en* +6. *hertz/en* If the *damping* keyword is not specified, the *viscoelastic* model is used by default. @@ -250,20 +250,20 @@ 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. -*enhooke* and *enhertz* models are useful for cases where 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) `, *enhooke* calculates the damping coefficient for the *hooke* model as: +*hooke/en* and *hertz/en* models are useful for cases where 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) `, *hooke/en* calculates the damping coefficient for the *hooke* model as: .. math:: \eta_n = \sqrt{\frac{4m_{eff}k_n}{1+\left( \frac{\pi}{\log(e)}\right)^2}}, -*enhertz* calculates the damping coefficient for the *hertz* and *hertz/material* models using: +*hertz/en* calculates the damping coefficient for the *hertz* and *hertz/material* models using: .. 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 these models calculate the damping coefficients by accounting for the effective mass, effective radius and pairwise overlaps (for *enhertz*), they accurately reproduce the specified coefficient of restitution for both monodisperse and polydisperse particle pairs. +where :math:`k_n = \frac{4}{3} E_{eff}` for the *hertz/material* model. Since these models calculate the damping coefficients by accounting for the effective mass, effective radius and pairwise overlaps (for *hertz/en*), they accurately reproduce 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/src/GRANULAR/gran_sub_mod_damping.cpp b/src/GRANULAR/gran_sub_mod_damping.cpp index e57972db767..12102575c39 100644 --- a/src/GRANULAR/gran_sub_mod_damping.cpp +++ b/src/GRANULAR/gran_sub_mod_damping.cpp @@ -141,44 +141,44 @@ double GranSubModDampingTsuji::calculate_forces() } /* ---------------------------------------------------------------------- - enhooke damping + hookeen damping ------------------------------------------------------------------------- */ -GranSubModDampingEnHooke::GranSubModDampingEnHooke(GranularModel *gm, LAMMPS *lmp) : +GranSubModDampingHookeEn::GranSubModDampingHookeEn(GranularModel *gm, LAMMPS *lmp) : GranSubModDamping(gm, lmp) { } -void GranSubModDampingEnHooke::init() +void GranSubModDampingHookeEn::init() { double cor = gm->normal_model->get_damp(); double logcor = log(cor); damp = -2*logcor/sqrt(PISQ + logcor*logcor); } -double GranSubModDampingEnHooke::calculate_forces() +double GranSubModDampingHookeEn::calculate_forces() { damp_prefactor = damp * sqrt(gm->meff * gm->Fnormal / gm->delta); return -damp_prefactor * gm->vnnr; } /* ---------------------------------------------------------------------- - enhertz damping + hertzen damping ------------------------------------------------------------------------- */ -GranSubModDampingEnHertz::GranSubModDampingEnHertz(GranularModel *gm, LAMMPS *lmp) : +GranSubModDampingHertzEn::GranSubModDampingHertzEn(GranularModel *gm, LAMMPS *lmp) : GranSubModDamping(gm, lmp) { } -void GranSubModDampingEnHertz::init() +void GranSubModDampingHertzEn::init() { double cor = gm->normal_model->get_damp(); double logcor = log(cor); damp = -ROOTTHREEBYTWO*TWOROOTFIVEBYSIX*logcor/sqrt(PISQ + logcor*logcor); } -double GranSubModDampingEnHertz::calculate_forces() +double GranSubModDampingHertzEn::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 33e88718fd3..9f37e144509 100644 --- a/src/GRANULAR/gran_sub_mod_damping.h +++ b/src/GRANULAR/gran_sub_mod_damping.h @@ -18,8 +18,8 @@ GranSubModStyle(velocity,GranSubModDampingVelocity,DAMPING); GranSubModStyle(mass_velocity,GranSubModDampingMassVelocity,DAMPING); GranSubModStyle(viscoelastic,GranSubModDampingViscoelastic,DAMPING); GranSubModStyle(tsuji,GranSubModDampingTsuji,DAMPING); -GranSubModStyle(enhooke,GranSubModDampingEnHooke,DAMPING); -GranSubModStyle(enhertz,GranSubModDampingEnHertz,DAMPING); +GranSubModStyle(hooke/en,GranSubModDampingHookeEn,DAMPING); +GranSubModStyle(hertz/en,GranSubModDampingHertzEn,DAMPING); // clang-format on #else @@ -88,18 +88,18 @@ namespace Granular_NS { /* ---------------------------------------------------------------------- */ - class GranSubModDampingEnHooke : public GranSubModDamping { + class GranSubModDampingHookeEn : public GranSubModDamping { public: - GranSubModDampingEnHooke(class GranularModel *, class LAMMPS *); + GranSubModDampingHookeEn(class GranularModel *, class LAMMPS *); void init() override; double calculate_forces() override; }; /* ---------------------------------------------------------------------- */ - class GranSubModDampingEnHertz : public GranSubModDamping { + class GranSubModDampingHertzEn : public GranSubModDamping { public: - GranSubModDampingEnHertz(class GranularModel *, class LAMMPS *); + GranSubModDampingHertzEn(class GranularModel *, class LAMMPS *); void init() override; double calculate_forces() override; }; diff --git a/src/GRANULAR/granular_model.cpp b/src/GRANULAR/granular_model.cpp index c1ad692fb3f..ed3c00866ec 100644 --- a/src/GRANULAR/granular_model.cpp +++ b/src/GRANULAR/granular_model.cpp @@ -238,6 +238,10 @@ void GranularModel::init() if (damping_model->name == "none") error->all(FLERR, "Must specify damping granular model"); if (tangential_model->name == "none") error->all(FLERR, "Must specify tangential granular model"); + //Check if correct damping model is being used with the normal model + if(normal_model->name =="hooke" && damping_model->name == "hertz/en") error->all(FLERR, "hooke should not be used with hertz/en damping model, please use hooke/en"); + if((normal_model->name =="hertz" || normal_model->name =="hertz/material") && damping_model->name == "hooke/en") error->all(FLERR, "hertz/material or hertz should not be used with hooke/en damping model, please use hertz/en"); + // Twisting, rolling, and heat are optional twisting_defined = rolling_defined = heat_defined = 1; if (twisting_model->name == "none") twisting_defined = 0; From 0344b6af707c7d36c199c3bfbb4546915c36627f Mon Sep 17 00:00:00 2001 From: Dhairya Date: Fri, 5 Apr 2024 10:14:07 -0500 Subject: [PATCH 08/13] updated the associated example file --- examples/granular/en_example/start.lammps | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/examples/granular/en_example/start.lammps b/examples/granular/en_example/start.lammps index 2bf0f54ddc0..7a2a56a0e91 100644 --- a/examples/granular/en_example/start.lammps +++ b/examples/granular/en_example/start.lammps @@ -9,9 +9,10 @@ read_data particles.dat 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 enhertz -# pair_coeff * * hertz 1e6 0.3 tangential mindlin 1e4 1.0 0.5 damping enhertz -# pair_coeff * * hooke 1e6 0.5 tangential mindlin 1 1.0 0.0 damping enhooke +pair_coeff * * hertz/material 1e6 0.8 0.4 tangential mindlin NULL 0.0 0.5 damping hertz/en +# pair_coeff * * hertz/material 1e6 0.8 0.4 tangential mindlin NULL 0.0 0.5 damping hooke/en #Should throw and error +# pair_coeff * * hertz 1e6 0.3 tangential mindlin 1e4 1.0 0.5 damping hertz/en +# pair_coeff * * hooke 1e6 0.5 tangential mindlin 1 1.0 0.0 damping hooke/en comm_modify vel yes timestep 1e-9 From 3d4bb44911d1440bda0472b32dbc4f4b75427f75 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Mon, 6 May 2024 16:39:29 -0600 Subject: [PATCH 09/13] Minor rearrangements to CoR, fix bug in granular single --- doc/src/pair_granular.rst | 21 +++-- .../data.particles} | 0 .../in.coeff_restitution} | 14 ++-- .../coeff_restitution/log.11Nov23.cor.g++.1 | 79 +++++++++++++++++++ src/GRANULAR/gran_sub_mod_damping.cpp | 45 ++++------- src/GRANULAR/gran_sub_mod_damping.h | 16 +--- src/GRANULAR/gran_sub_mod_normal.h | 3 +- src/GRANULAR/granular_model.cpp | 4 - src/GRANULAR/pair_granular.cpp | 2 +- 9 files changed, 120 insertions(+), 64 deletions(-) rename examples/granular/{en_example/particles.dat => coeff_restitution/data.particles} (100%) rename examples/granular/{en_example/start.lammps => coeff_restitution/in.coeff_restitution} (56%) create mode 100644 examples/granular/coeff_restitution/log.11Nov23.cor.g++.1 diff --git a/doc/src/pair_granular.rst b/doc/src/pair_granular.rst index 9eb34c9de84..597869a05b4 100644 --- a/doc/src/pair_granular.rst +++ b/doc/src/pair_granular.rst @@ -187,8 +187,7 @@ for the damping model currently supported are: 2. *mass_velocity* 3. *viscoelastic* 4. *tsuji* -5. *hooke/en* -6. *hertz/en* +5. *coeff_restitution* If the *damping* keyword is not specified, the *viscoelastic* model is used by default. @@ -250,20 +249,28 @@ 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. -*hooke/en* and *hertz/en* models are useful for cases where 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) `, *hooke/en* calculates the damping coefficient for the *hooke* model as: +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}}, + \eta_n = \sqrt{\frac{4m_{eff}k_n}{1+\left( \frac{\pi}{\log(e)}\right)^2}} , -*hertz/en* calculates the damping coefficient for the *hertz* and *hertz/material* models using: +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 these models calculate the damping coefficients by accounting for the effective mass, effective radius and pairwise overlaps (for *hertz/en*), they accurately reproduce the specified coefficient of restitution for both monodisperse and polydisperse particle pairs. +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/examples/granular/en_example/particles.dat b/examples/granular/coeff_restitution/data.particles similarity index 100% rename from examples/granular/en_example/particles.dat rename to examples/granular/coeff_restitution/data.particles diff --git a/examples/granular/en_example/start.lammps b/examples/granular/coeff_restitution/in.coeff_restitution similarity index 56% rename from examples/granular/en_example/start.lammps rename to examples/granular/coeff_restitution/in.coeff_restitution index 7a2a56a0e91..e441ed67a75 100644 --- a/examples/granular/en_example/start.lammps +++ b/examples/granular/coeff_restitution/in.coeff_restitution @@ -1,26 +1,24 @@ -units si +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 particles.dat add append +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 hertz/en -# pair_coeff * * hertz/material 1e6 0.8 0.4 tangential mindlin NULL 0.0 0.5 damping hooke/en #Should throw and error -# pair_coeff * * hertz 1e6 0.3 tangential mindlin 1e4 1.0 0.5 damping hertz/en -# pair_coeff * * hooke 1e6 0.5 tangential mindlin 1 1.0 0.0 damping hooke/en +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 +#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/coeff_restitution/log.11Nov23.cor.g++.1 b/examples/granular/coeff_restitution/log.11Nov23.cor.g++.1 new file mode 100644 index 00000000000..b4f53c84da8 --- /dev/null +++ b/examples/granular/coeff_restitution/log.11Nov23.cor.g++.1 @@ -0,0 +1,79 @@ +LAMMPS (21 Nov 2023 - Development - ) +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.06865 on 1 procs for 10000000 steps with 2 atoms + +99.7% CPU use with 1 MPI tasks x no OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 0.55389 | 0.55389 | 0.55389 | 0.0 | 10.93 +Neigh | 0.00022182 | 0.00022182 | 0.00022182 | 0.0 | 0.00 +Comm | 1.9988 | 1.9988 | 1.9988 | 0.0 | 39.43 +Output | 0.00023837 | 0.00023837 | 0.00023837 | 0.0 | 0.00 +Modify | 1.26 | 1.26 | 1.26 | 0.0 | 24.86 +Other | | 1.255 | | | 24.77 + +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:05 diff --git a/src/GRANULAR/gran_sub_mod_damping.cpp b/src/GRANULAR/gran_sub_mod_damping.cpp index 12102575c39..be87c70cd48 100644 --- a/src/GRANULAR/gran_sub_mod_damping.cpp +++ b/src/GRANULAR/gran_sub_mod_damping.cpp @@ -16,19 +16,20 @@ #include "gran_sub_mod_normal.h" #include "granular_model.h" #include "math_special.h" +#include "math_const.h" #include using namespace LAMMPS_NS; using namespace Granular_NS; +using namespace MathConst; using MathSpecial::cube; using MathSpecial::powint; using MathSpecial::square; -static constexpr double PISQ = 9.86960440108935799230; -static constexpr double TWOROOTFIVEBYSIX = 1.82574185835055380345; -static constexpr double ROOTTHREEBYTWO = 1.22474487139158894067; +static constexpr double TWOROOTFIVEBYSIX = 1.82574185835055380345; // 2/sqrt(5/6) +static constexpr double ROOTTHREEBYTWO = 1.22474487139158894067; // sqrt(3/2) /* ---------------------------------------------------------------------- Default damping model @@ -141,44 +142,28 @@ double GranSubModDampingTsuji::calculate_forces() } /* ---------------------------------------------------------------------- - hookeen damping + Coefficient of restitution damping ------------------------------------------------------------------------- */ -GranSubModDampingHookeEn::GranSubModDampingHookeEn(GranularModel *gm, LAMMPS *lmp) : +GranSubModDampingCoeffRestitution::GranSubModDampingCoeffRestitution(GranularModel *gm, LAMMPS *lmp) : GranSubModDamping(gm, lmp) { } -void GranSubModDampingHookeEn::init() +void GranSubModDampingCoeffRestitution::init() { + // Calculate prefactor, assume Hertzian as default double cor = gm->normal_model->get_damp(); double logcor = log(cor); - damp = -2*logcor/sqrt(PISQ + logcor*logcor); + 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 GranSubModDampingHookeEn::calculate_forces() -{ - damp_prefactor = damp * sqrt(gm->meff * gm->Fnormal / gm->delta); - return -damp_prefactor * gm->vnnr; -} - -/* ---------------------------------------------------------------------- - hertzen damping -------------------------------------------------------------------------- */ - -GranSubModDampingHertzEn::GranSubModDampingHertzEn(GranularModel *gm, LAMMPS *lmp) : - GranSubModDamping(gm, lmp) -{ -} - -void GranSubModDampingHertzEn::init() -{ - double cor = gm->normal_model->get_damp(); - double logcor = log(cor); - damp = -ROOTTHREEBYTWO*TWOROOTFIVEBYSIX*logcor/sqrt(PISQ + logcor*logcor); -} - -double GranSubModDampingHertzEn::calculate_forces() +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 9f37e144509..c931e385cca 100644 --- a/src/GRANULAR/gran_sub_mod_damping.h +++ b/src/GRANULAR/gran_sub_mod_damping.h @@ -18,8 +18,7 @@ GranSubModStyle(velocity,GranSubModDampingVelocity,DAMPING); GranSubModStyle(mass_velocity,GranSubModDampingMassVelocity,DAMPING); GranSubModStyle(viscoelastic,GranSubModDampingViscoelastic,DAMPING); GranSubModStyle(tsuji,GranSubModDampingTsuji,DAMPING); -GranSubModStyle(hooke/en,GranSubModDampingHookeEn,DAMPING); -GranSubModStyle(hertz/en,GranSubModDampingHertzEn,DAMPING); +GranSubModStyle(coeff_restitution,GranSubModDampingCoeffRestitution,DAMPING); // clang-format on #else @@ -88,18 +87,9 @@ namespace Granular_NS { /* ---------------------------------------------------------------------- */ - class GranSubModDampingHookeEn : public GranSubModDamping { + class GranSubModDampingCoeffRestitution : public GranSubModDamping { public: - GranSubModDampingHookeEn(class GranularModel *, class LAMMPS *); - void init() override; - double calculate_forces() override; - }; - - /* ---------------------------------------------------------------------- */ - - class GranSubModDampingHertzEn : public GranSubModDamping { - public: - GranSubModDampingHertzEn(class GranularModel *, class LAMMPS *); + GranSubModDampingCoeffRestitution(class GranularModel *, class LAMMPS *); void init() override; double calculate_forces() override; }; 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/granular_model.cpp b/src/GRANULAR/granular_model.cpp index ed3c00866ec..c1ad692fb3f 100644 --- a/src/GRANULAR/granular_model.cpp +++ b/src/GRANULAR/granular_model.cpp @@ -238,10 +238,6 @@ void GranularModel::init() if (damping_model->name == "none") error->all(FLERR, "Must specify damping granular model"); if (tangential_model->name == "none") error->all(FLERR, "Must specify tangential granular model"); - //Check if correct damping model is being used with the normal model - if(normal_model->name =="hooke" && damping_model->name == "hertz/en") error->all(FLERR, "hooke should not be used with hertz/en damping model, please use hooke/en"); - if((normal_model->name =="hertz" || normal_model->name =="hertz/material") && damping_model->name == "hooke/en") error->all(FLERR, "hertz/material or hertz should not be used with hooke/en damping model, please use hertz/en"); - // Twisting, rolling, and heat are optional twisting_defined = rolling_defined = heat_defined = 1; if (twisting_model->name == "none") twisting_defined = 0; diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index 119feb1c380..480e9084875 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -770,7 +770,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]; From 33525de5985434bdf8fed7801fefaaabba2ef2b2 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Wed, 8 May 2024 16:39:25 -0600 Subject: [PATCH 10/13] fix incomplete header info on command syntax --- doc/src/write_data.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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 From d4d4c4857466fa60d13998a6ed3db23bcfe6ebee Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 13 May 2024 13:43:50 -0400 Subject: [PATCH 11/13] spelling fixes --- doc/src/replicate.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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. From 4e7bddaa0b3286658dcbaa1c062b2b79ba38ae4e Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 13 May 2024 13:49:59 -0400 Subject: [PATCH 12/13] remove unused variables --- src/EXTRA-FIX/fix_deform_pressure.cpp | 2 -- src/ML-IAP/mliap_unified.cpp | 2 -- src/create_atoms.cpp | 2 +- 3 files changed, 1 insertion(+), 5 deletions(-) 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/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); From 74f29ba2f3fff1b5f3af958b52915b4a9127be12 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 13 May 2024 14:07:24 -0400 Subject: [PATCH 13/13] update example --- .../{coeff_restitution => }/data.particles | 0 .../in.coeff_restitution => in.restitution} | 0 ...or.g++.1 => log.13May23.restitution.g++.1} | 25 ++++++++++--------- 3 files changed, 13 insertions(+), 12 deletions(-) rename examples/granular/{coeff_restitution => }/data.particles (100%) rename examples/granular/{coeff_restitution/in.coeff_restitution => in.restitution} (100%) rename examples/granular/{coeff_restitution/log.11Nov23.cor.g++.1 => log.13May23.restitution.g++.1} (74%) diff --git a/examples/granular/coeff_restitution/data.particles b/examples/granular/data.particles similarity index 100% rename from examples/granular/coeff_restitution/data.particles rename to examples/granular/data.particles diff --git a/examples/granular/coeff_restitution/in.coeff_restitution b/examples/granular/in.restitution similarity index 100% rename from examples/granular/coeff_restitution/in.coeff_restitution rename to examples/granular/in.restitution diff --git a/examples/granular/coeff_restitution/log.11Nov23.cor.g++.1 b/examples/granular/log.13May23.restitution.g++.1 similarity index 74% rename from examples/granular/coeff_restitution/log.11Nov23.cor.g++.1 rename to examples/granular/log.13May23.restitution.g++.1 index b4f53c84da8..e51709d10dc 100644 --- a/examples/granular/coeff_restitution/log.11Nov23.cor.g++.1 +++ b/examples/granular/log.13May23.restitution.g++.1 @@ -1,4 +1,5 @@ -LAMMPS (21 Nov 2023 - Development - ) +LAMMPS (17 Apr 2024 - Development - patch_17Apr2024-93-g4e7bddaa0b) + using 1 OpenMP thread(s) per MPI task units si atom_style sphere @@ -29,8 +30,8 @@ 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 +#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 @@ -51,19 +52,19 @@ 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.06865 on 1 procs for 10000000 steps with 2 atoms +Loop time of 5.99782 on 1 procs for 10000000 steps with 2 atoms -99.7% CPU use with 1 MPI tasks x no OpenMP threads +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.55389 | 0.55389 | 0.55389 | 0.0 | 10.93 -Neigh | 0.00022182 | 0.00022182 | 0.00022182 | 0.0 | 0.00 -Comm | 1.9988 | 1.9988 | 1.9988 | 0.0 | 39.43 -Output | 0.00023837 | 0.00023837 | 0.00023837 | 0.0 | 0.00 -Modify | 1.26 | 1.26 | 1.26 | 0.0 | 24.86 -Other | | 1.255 | | | 24.77 +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 @@ -76,4 +77,4 @@ Total # of neighbors = 0 Ave neighs/atom = 0 Neighbor list builds = 14 Dangerous builds = 0 -Total wall time: 0:00:05 +Total wall time: 0:00:06