diff --git a/doc/src/bond_bpm_rotational.rst b/doc/src/bond_bpm_rotational.rst index ba93d679bab..ca12d86cccb 100644 --- a/doc/src/bond_bpm_rotational.rst +++ b/doc/src/bond_bpm_rotational.rst @@ -24,14 +24,17 @@ Syntax *x, y, z* = the center of mass position of the 2 atoms when the bond broke (distance units) *x/ref, y/ref, z/ref* = the initial center of mass position of the 2 atoms (distance units) - *overlay/pair* value = none + *overlay/pair* value = *yes* or *no* bonded particles will still interact with pair forces *smooth* value = *yes* or *no* smooths bond forces near the breaking point - *break/no* - indicates that bonds should not break during a run + *normalize* value = *yes* or *no* + normalizes normal and shear forces by the reference length + + *break* value = *yes* or *no* + indicates whether bonds break during a run Examples """""""" @@ -136,16 +139,19 @@ or :doc:`read_restart ` commands: * :math:`\gamma_r` (force*distance/velocity units) * :math:`\gamma_t` (force*distance/velocity units) +However, the *normalize* option will normalize the radial and shear forces +by :math:`r_0` such that :math:`k_r` and :math:`k_s` are unit less. + By default, pair forces are not calculated between bonded particles. Pair forces can alternatively be overlaid on top of bond forces using -the *overlay/pair* keyword. These settings require specific +the *overlay/pair* option. These settings require specific :doc:`special_bonds ` settings described in the restrictions. Further details can be found in the `:doc: how to ` page on BPMs. .. versionadded:: 28Mar2023 -If the *break/no* keyword is used, then LAMMPS assumes bonds should not break +If the *break* option is used, then LAMMPS assumes bonds should not break during a simulation run. This will prevent some unnecessary calculation. However, if a bond does break, it will trigger an error. @@ -251,7 +257,7 @@ Related commands Default """"""" -The option defaults are *smooth* = *yes* +The option defaults are *overlay/pair* = *no*, *smooth* = *yes*, *normalize* = *no*, and *break* = *yes* ---------- diff --git a/doc/src/bond_bpm_spring.rst b/doc/src/bond_bpm_spring.rst index 5762dbe2089..d89035dcad5 100644 --- a/doc/src/bond_bpm_spring.rst +++ b/doc/src/bond_bpm_spring.rst @@ -24,14 +24,17 @@ Syntax *x, y, z* = the center of mass position of the 2 atoms when the bond broke (distance units) *x/ref, y/ref, z/ref* = the initial center of mass position of the 2 atoms (distance units) - *overlay/pair* value = none + *overlay/pair* value = *yes* or *no* bonded particles will still interact with pair forces *smooth* value = *yes* or *no* smooths bond forces near the breaking point - *break/no* - indicates that bonds should not break during a run + *normalize* value = *yes* or *no* + normalizes bond forces by the reference length + + *break* value = *yes* or *no* + indicates whether bonds break during a run Examples """""""" @@ -66,7 +69,7 @@ particles based on a model described by Clemmer and Robbins F = k (r - r_0) w -where :math:`k_r` is a stiffness, :math:`r` is the current distance +where :math:`k` is a stiffness, :math:`r` is the current distance and :math:`r_0` is the initial distance between the two particles, and :math:`w` is an optional smoothing factor discussed below. Bonds will break at a strain of :math:`\epsilon_c`. This is done by setting by @@ -102,16 +105,19 @@ the data file or restart files read by the :doc:`read_data * :math:`\epsilon_c` (unit less) * :math:`\gamma` (force/velocity units) +However, the *normalize* option will normalize the elastic bond force by +:math:`r_0` such that :math:`k` is unit less. + By default, pair forces are not calculated between bonded particles. Pair forces can alternatively be overlaid on top of bond forces using -the *overlay/pair* keyword. These settings require specific +the *overlay/pair* option. These settings require specific :doc:`special_bonds ` settings described in the restrictions. Further details can be found in the `:doc: how to ` page on BPMs. .. versionadded:: 28Mar2023 -If the *break/no* keyword is used, then LAMMPS assumes bonds should not break +If the *break* option is used, then LAMMPS assumes bonds should not break during a simulation run. This will prevent some unnecessary calculation. However, if a bond does break, it will trigger an error. @@ -206,7 +212,7 @@ Related commands Default """"""" -The option defaults are *smooth* = *yes* +The option defaults are *overlay/pair* = *no*, *smooth* = *yes*, *normalize* = *no*, and *break* = *yes* ---------- diff --git a/doc/src/compute_bond_local.rst b/doc/src/compute_bond_local.rst index f3fb752ebe7..10e86bbe442 100644 --- a/doc/src/compute_bond_local.rst +++ b/doc/src/compute_bond_local.rst @@ -76,7 +76,10 @@ The value *force* is the magnitude of the force acting between the pair of atoms in the bond. The values *fx*, *fy*, and *fz* are the xyz components of -*force* between the pair of atoms in the bond. +*force* between the pair of atoms in the bond. For bond styles that apply +non-central forces, such as :doc:`bond_style bpm/rotational +`, these values only include the :math:`(x,y,z)` +components of the normal force component. The remaining properties are all computed for motion of the two atoms relative to the center of mass (COM) velocity of the 2 atoms in the diff --git a/doc/src/compute_fabric.rst b/doc/src/compute_fabric.rst index f1a5d3d7f9f..b38ffafa48a 100644 --- a/doc/src/compute_fabric.rst +++ b/doc/src/compute_fabric.rst @@ -146,13 +146,13 @@ m to :math:`M` (inclusive). A middle asterisk means all types from m to n Output info """"""""""" -This compute calculates a local vector of doubles and a scalar. The vector -stores the unique components of the first requested tensor in the order -:math:`xx`, :math:`yy`, :math:`zz`, :math:`xy`, :math:`xz`, :math:`yz` -followed by the same components for all subsequent tensors. +This compute calculates a global vector of doubles and a global scalar. The +vector stores the unique components of the first requested tensor in the +order :math:`xx`, :math:`yy`, :math:`zz`, :math:`xy`, :math:`xz`, +:math:`yz` followed by the same components for all subsequent tensors. The length of the vector is therefore six times the number of requested -tensors. The scalar output is the number of pairwise interactions included in -the calculation of the fabric tensor. +tensors. The scalar output is the number of pairwise interactions included +in the calculation of the fabric tensor. Restrictions """""""""""" diff --git a/doc/src/compute_pair_local.rst b/doc/src/compute_pair_local.rst index dace280dee1..31209f63f49 100644 --- a/doc/src/compute_pair_local.rst +++ b/doc/src/compute_pair_local.rst @@ -66,7 +66,9 @@ The value *eng* is the interaction energy for the pair of atoms. The value *force* is the force acting between the pair of atoms, which is positive for a repulsive force and negative for an attractive force. The values *fx*, *fy*, and *fz* are the :math:`(x,y,z)` components of -*force* on atom I. +*force* on atom I. For pair styles that apply non-central forces, +such as :doc:`granular pair styles `, these values only include +the :math:`(x,y,z)` components of the normal force component. A pair style may define additional pairwise quantities which can be accessed as *p1* to *pN*, where :math:`N` is defined by the pair style. diff --git a/doc/src/pair_granular.rst b/doc/src/pair_granular.rst index 21afc1b4fc7..37ff8fcbfec 100644 --- a/doc/src/pair_granular.rst +++ b/doc/src/pair_granular.rst @@ -736,7 +736,7 @@ or .. math:: - E_{eff,ij} = \frac{E_{ij}}{2(1-\nu_{ij})} + E_{eff,ij} = \frac{E_{ij}}{2(1-\nu_{ij}^2)} These pair styles write their information to :doc:`binary restart files `, so a pair_style command does not need to be specified in an input script that reads a restart file. diff --git a/src/BPM/bond_bpm.cpp b/src/BPM/bond_bpm.cpp index a68aea3e683..da189ccaf54 100644 --- a/src/BPM/bond_bpm.cpp +++ b/src/BPM/bond_bpm.cpp @@ -187,10 +187,12 @@ void BondBPM::settings(int narg, char **arg) iarg++; } } else if (strcmp(arg[iarg], "overlay/pair") == 0) { - overlay_flag = 1; + if (iarg + 1 > narg) error->all(FLERR, "Illegal bond bpm command, missing option for overlay/pair"); + overlay_flag = utils::logical(FLERR, arg[iarg + 1], false, lmp); iarg++; - } else if (strcmp(arg[iarg], "break/no") == 0) { - break_flag = 0; + } else if (strcmp(arg[iarg], "break") == 0) { + if (iarg + 1 > narg) error->all(FLERR, "Illegal bond bpm command, missing option for break"); + break_flag = utils::logical(FLERR, arg[iarg + 1], false, lmp); iarg++; } else { leftover_iarg.push_back(iarg); diff --git a/src/BPM/bond_bpm_rotational.cpp b/src/BPM/bond_bpm_rotational.cpp index ac29c0e376f..ffb0d9521d0 100644 --- a/src/BPM/bond_bpm_rotational.cpp +++ b/src/BPM/bond_bpm_rotational.cpp @@ -50,6 +50,7 @@ BondBPMRotational::BondBPMRotational(LAMMPS *_lmp) : { partial_flag = 1; smooth_flag = 1; + normalize_flag = 0; single_extra = 7; svector = new double[7]; @@ -203,6 +204,13 @@ double BondBPMRotational::elastic_forces(int i1, int i2, int type, double r_mag, double Ts[3], Tb[3], Tt[3], Tbp[3], Ttp[3], Tsp[3], T_rot[3], Ttmp[3]; double **quat = atom->quat; + double r0_mag_inv = 1.0 / r0_mag; + double Kr_type = Kr[type]; + double Ks_type = Ks[type]; + if (normalize_flag) { + Kr_type *= r0_mag_inv; + Ks_type *= r0_mag_inv; + } q1[0] = quat[i1][0]; q1[1] = quat[i1][1]; @@ -217,26 +225,26 @@ double BondBPMRotational::elastic_forces(int i1, int i2, int type, double r_mag, // Calculate normal forces, rb = bond vector in particle 1's frame MathExtra::qconjugate(q2, q2inv); MathExtra::quatrotvec(q2inv, r, rb); - Fr = Kr[type] * (r_mag - r0_mag); + Fr = Kr_type * (r_mag - r0_mag); MathExtra::scale3(Fr * r_mag_inv, rb, F_rot); // Calculate forces due to tangential displacements (no rotation) r0_dot_rb = MathExtra::dot3(r0, rb); - c = r0_dot_rb * r_mag_inv / r0_mag; + c = r0_dot_rb * r_mag_inv * r0_mag_inv; gamma = acos_limit(c); MathExtra::cross3(rb, r0, rb_x_r0); MathExtra::cross3(rb, rb_x_r0, s); MathExtra::norm3(s); - MathExtra::scale3(Ks[type] * r_mag * gamma, s, Fs); + MathExtra::scale3(Ks_type * r_mag * gamma, s, Fs); // Calculate torque due to tangential displacements MathExtra::cross3(r0, rb, t); MathExtra::norm3(t); - MathExtra::scale3(0.5 * r_mag * Ks[type] * r_mag * gamma, t, Ts); + MathExtra::scale3(0.5 * r_mag * Ks_type * r_mag * gamma, t, Ts); // Relative rotation force/torque // Use representation of X'Y'Z' rotations from Wang, Mora 2009 @@ -316,12 +324,12 @@ double BondBPMRotational::elastic_forces(int i1, int i2, int type, double r_mag, Ttp[1] = 0.0; Ttp[2] = Kt[type] * psi; - Fsp[0] = -0.5 * Ks[type] * r_mag * theta * cos_phi; - Fsp[1] = -0.5 * Ks[type] * r_mag * theta * sin_phi; + Fsp[0] = -0.5 * Ks_type * r_mag * theta * cos_phi; + Fsp[1] = -0.5 * Ks_type * r_mag * theta * sin_phi; Fsp[2] = 0.0; - Tsp[0] = 0.25 * Ks[type] * r_mag * r_mag * theta * sin_phi; - Tsp[1] = -0.25 * Ks[type] * r_mag * r_mag * theta * cos_phi; + Tsp[0] = 0.25 * Ks_type * r_mag * r_mag * theta * sin_phi; + Tsp[1] = -0.25 * Ks_type * r_mag * r_mag * theta * cos_phi; Tsp[2] = 0.0; // Rotate forces/torques back to 1st particle's frame @@ -667,6 +675,10 @@ void BondBPMRotational::settings(int narg, char **arg) if (iarg + 1 > narg) error->all(FLERR, "Illegal bond bpm command, missing option for smooth"); smooth_flag = utils::logical(FLERR, arg[iarg + 1], false, lmp); i += 1; + } else if (strcmp(arg[iarg], "normalize") == 0) { + if (iarg + 1 > narg) error->all(FLERR, "Illegal bond bpm command, missing option for normalize"); + normalize_flag = utils::logical(FLERR, arg[iarg + 1], false, lmp); + i += 1; } else { error->all(FLERR, "Illegal bond bpm command, invalid argument {}", arg[iarg]); } @@ -793,7 +805,7 @@ double BondBPMRotational::single(int type, double rsq, int i, int j, double &ffo double breaking = elastic_forces(i, j, type, r_mag, r0_mag, r_mag_inv, rhat, r, r0, force1on2, torque1on2, torque2on1); damping_forces(i, j, type, rhat, r, force1on2, torque1on2, torque2on1); - fforce = MathExtra::dot3(force1on2, r); + fforce = MathExtra::dot3(force1on2, rhat); fforce *= -1; double smooth = 1.0; diff --git a/src/BPM/bond_bpm_rotational.h b/src/BPM/bond_bpm_rotational.h index 0fb38e73437..9fdc418d2d2 100644 --- a/src/BPM/bond_bpm_rotational.h +++ b/src/BPM/bond_bpm_rotational.h @@ -41,7 +41,7 @@ class BondBPMRotational : public BondBPM { protected: double *Kr, *Ks, *Kt, *Kb, *gnorm, *gslide, *groll, *gtwist; double *Fcr, *Fcs, *Tct, *Tcb; - int smooth_flag; + int smooth_flag, normalize_flag; double elastic_forces(int, int, int, double, double, double, double *, double *, double *, double *, double *, double *); diff --git a/src/BPM/bond_bpm_spring.cpp b/src/BPM/bond_bpm_spring.cpp index ed4e71daf1b..37b79f93fba 100644 --- a/src/BPM/bond_bpm_spring.cpp +++ b/src/BPM/bond_bpm_spring.cpp @@ -37,6 +37,7 @@ BondBPMSpring::BondBPMSpring(LAMMPS *_lmp) : { partial_flag = 1; smooth_flag = 1; + normalize_flag = 0; single_extra = 1; svector = new double[1]; @@ -190,7 +191,10 @@ void BondBPMSpring::compute(int eflag, int vflag) } rinv = 1.0 / r; - fbond = k[type] * (r0 - r); + if (normalize_flag) + fbond = -k[type] * e; + else + fbond = k[type] * (r0 - r); delvx = v[i1][0] - v[i2][0]; delvy = v[i1][1] - v[i2][1]; @@ -302,6 +306,10 @@ void BondBPMSpring::settings(int narg, char **arg) if (iarg + 1 > narg) error->all(FLERR, "Illegal bond bpm command, missing option for smooth"); smooth_flag = utils::logical(FLERR, arg[iarg + 1], false, lmp); i += 1; + } else if (strcmp(arg[iarg], "normalize") == 0) { + if (iarg + 1 > narg) error->all(FLERR, "Illegal bond bpm command, missing option for normalize"); + normalize_flag = utils::logical(FLERR, arg[iarg + 1], false, lmp); + i += 1; } else { error->all(FLERR, "Illegal bond bpm command, invalid argument {}", arg[iarg]); } @@ -376,7 +384,11 @@ double BondBPMSpring::single(int type, double rsq, int i, int j, double &fforce) double r = sqrt(rsq); double rinv = 1.0 / r; - fforce = k[type] * (r0 - r); + + if (normalize_flag) + fforce = k[type] * (r0 - r) / r0; + else + fforce = k[type] * (r0 - r); double **x = atom->x; double **v = atom->v; diff --git a/src/BPM/bond_bpm_spring.h b/src/BPM/bond_bpm_spring.h index 409469bef3e..93f4b49a26a 100644 --- a/src/BPM/bond_bpm_spring.h +++ b/src/BPM/bond_bpm_spring.h @@ -40,7 +40,7 @@ class BondBPMSpring : public BondBPM { protected: double *k, *ecrit, *gamma; - int smooth_flag; + int smooth_flag, normalize_flag; void allocate(); void store_data(); diff --git a/src/GRANULAR/gran_sub_mod_normal.cpp b/src/GRANULAR/gran_sub_mod_normal.cpp index fcd9cdd27fb..ffc18b8c327 100644 --- a/src/GRANULAR/gran_sub_mod_normal.cpp +++ b/src/GRANULAR/gran_sub_mod_normal.cpp @@ -150,6 +150,7 @@ GranSubModNormalHertzMaterial::GranSubModNormalHertzMaterial(GranularModel *gm, material_properties = 1; num_coeffs = 3; contact_radius_flag = 1; + mixed_coefficients = 0; } /* ---------------------------------------------------------------------- */ @@ -159,10 +160,12 @@ void GranSubModNormalHertzMaterial::coeffs_to_local() Emod = coeffs[0]; damp = coeffs[1]; poiss = coeffs[2]; - if (gm->contact_type == PAIR) { - k = FOURTHIRDS * mix_stiffnessE(Emod, Emod, poiss, poiss); - } else { - k = FOURTHIRDS * mix_stiffnessE_wall(Emod, poiss); + if (!mixed_coefficients) { + if (gm->contact_type == PAIR) { + k = FOURTHIRDS * mix_stiffnessE(Emod, Emod, poiss, poiss); + } else { + k = FOURTHIRDS * mix_stiffnessE_wall(Emod, poiss); + } } if (Emod < 0.0 || damp < 0.0) error->all(FLERR, "Illegal Hertz material normal model"); @@ -175,6 +178,10 @@ void GranSubModNormalHertzMaterial::mix_coeffs(double *icoeffs, double *jcoeffs) coeffs[0] = mix_stiffnessE(icoeffs[0], jcoeffs[0], icoeffs[2], jcoeffs[2]); coeffs[1] = mix_geom(icoeffs[1], jcoeffs[1]); coeffs[2] = mix_geom(icoeffs[2], jcoeffs[2]); + + k = FOURTHIRDS * coeffs[0]; + mixed_coefficients = 1; + coeffs_to_local(); } @@ -188,6 +195,7 @@ GranSubModNormalDMT::GranSubModNormalDMT(GranularModel *gm, LAMMPS *lmp) : GranS cohesive_flag = 1; num_coeffs = 4; contact_radius_flag = 1; + mixed_coefficients = 0; } /* ---------------------------------------------------------------------- */ @@ -198,10 +206,13 @@ void GranSubModNormalDMT::coeffs_to_local() damp = coeffs[1]; poiss = coeffs[2]; cohesion = coeffs[3]; - if (gm->contact_type == PAIR) { - k = FOURTHIRDS * mix_stiffnessE(Emod, Emod, poiss, poiss); - } else { - k = FOURTHIRDS * mix_stiffnessE_wall(Emod, poiss); + + if (!mixed_coefficients) { + if (gm->contact_type == PAIR) { + k = FOURTHIRDS * mix_stiffnessE(Emod, Emod, poiss, poiss); + } else { + k = FOURTHIRDS * mix_stiffnessE_wall(Emod, poiss); + } } if (Emod < 0.0 || damp < 0.0) error->all(FLERR, "Illegal DMT normal model"); @@ -215,6 +226,10 @@ void GranSubModNormalDMT::mix_coeffs(double *icoeffs, double *jcoeffs) coeffs[1] = mix_geom(icoeffs[1], jcoeffs[1]); coeffs[2] = mix_geom(icoeffs[2], jcoeffs[2]); coeffs[3] = mix_geom(icoeffs[3], jcoeffs[3]); + + k = FOURTHIRDS * coeffs[0]; + mixed_coefficients = 1; + coeffs_to_local(); } @@ -246,6 +261,7 @@ GranSubModNormalJKR::GranSubModNormalJKR(GranularModel *gm, LAMMPS *lmp) : GranS beyond_contact = 1; num_coeffs = 4; contact_radius_flag = 1; + mixed_coefficients = 0; } /* ---------------------------------------------------------------------- */ @@ -257,10 +273,12 @@ void GranSubModNormalJKR::coeffs_to_local() poiss = coeffs[2]; cohesion = coeffs[3]; - if (gm->contact_type == PAIR) { - Emix = mix_stiffnessE(Emod, Emod, poiss, poiss); - } else { - Emix = mix_stiffnessE_wall(Emod, poiss); + if (!mixed_coefficients) { + if (gm->contact_type == PAIR) { + Emix = mix_stiffnessE(Emod, Emod, poiss, poiss); + } else { + Emix = mix_stiffnessE_wall(Emod, poiss); + } } k = FOURTHIRDS * Emix; @@ -276,6 +294,10 @@ void GranSubModNormalJKR::mix_coeffs(double *icoeffs, double *jcoeffs) coeffs[1] = mix_geom(icoeffs[1], jcoeffs[1]); coeffs[2] = mix_geom(icoeffs[2], jcoeffs[2]); coeffs[3] = mix_geom(icoeffs[3], jcoeffs[3]); + + Emix = coeffs[0]; + mixed_coefficients = 1; + coeffs_to_local(); } diff --git a/src/GRANULAR/gran_sub_mod_normal.h b/src/GRANULAR/gran_sub_mod_normal.h index 7fb4a259e05..c1511fdfa5b 100644 --- a/src/GRANULAR/gran_sub_mod_normal.h +++ b/src/GRANULAR/gran_sub_mod_normal.h @@ -29,6 +29,7 @@ GranSubModStyle(jkr,GranSubModNormalJKR,NORMAL); namespace LAMMPS_NS { namespace Granular_NS { + class GranSubModNormal : public GranSubMod { public: GranSubModNormal(class GranularModel *, class LAMMPS *); @@ -92,6 +93,8 @@ namespace Granular_NS { GranSubModNormalHertzMaterial(class GranularModel *, class LAMMPS *); void coeffs_to_local() override; void mix_coeffs(double *, double *) override; + private: + int mixed_coefficients; }; /* ---------------------------------------------------------------------- */ @@ -107,6 +110,7 @@ namespace Granular_NS { protected: double k, cohesion; double F_pulloff, Fne; + int mixed_coefficients; }; /* ---------------------------------------------------------------------- */ @@ -125,6 +129,7 @@ namespace Granular_NS { protected: double k, cohesion; double Emix, F_pulloff, Fne; + int mixed_coefficients; }; } // namespace Granular_NS