diff --git a/doc/src/pair_granular.rst b/doc/src/pair_granular.rst index 9eb34c9de8..597869a05b 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 7a2a56a0e9..e441ed67a7 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 0000000000..b4f53c84da --- /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 12102575c3..be87c70cd4 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 9f37e14450..c931e385cc 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 6f2f3aabbe..c1ed36b6e1 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 ed3c00866e..c1ad692fb3 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 119feb1c38..480e908487 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];