apply changes suggested by Andrew Santos (747bfa31) · Commits · 郑智淋 / lammps

doc/src/pair_granular.txt

+7 −7

Original line number	Diff line number	Diff line
		@@ -177,7 +177,7 @@ following general form:
		\end\{equation\}

		Here, \(\mathbf\{v\}_\{n,rel\} = (\mathbf\{v\}_j - \mathbf\{v\}_i)
		\cdot \mathbf\{n\}\) is the component of relative velocity along
		\cdot \mathbf\{n\} \mathbf\{n\}\) is the component of relative velocity along
		\(\mathbf\{n\}\).

		The optional {damping} keyword to the {pair_coeff} command followed by
		@@ -299,8 +299,8 @@ the normal damping \(\eta_n\) (see above):
		\eta_t = -x_\{\gamma,t\} \eta_n
		\end\{equation\}

		The normal damping prefactor \(\eta_n\) is determined by the choice of
		the {damping} keyword, as discussed above. Thus, the {damping}
		The normal damping prefactor \(\eta_n\) is determined by the choice
		of the {damping} keyword, as discussed above. Thus, the {damping}
		keyword also affects the tangential damping. The parameter
		\(x_\{\gamma,t\}\) is a scaling coefficient. Several works in the
		literature use \(x_\{\gamma,t\} = 1\) ("Marshall"_#Marshall2009,
		@@ -308,9 +308,9 @@ literature use \(x_\{\gamma,t\} = 1\) ("Marshall"_#Marshall2009,
		tangential velocity at the point of contact is given by
		\(\mathbf\{v\}_\{t, rel\} = \mathbf\{v\}_\{t\} - (R_i\Omega_i +
		R_j\Omega_j) \times \mathbf\{n\}\), where \(\mathbf\{v\}_\{t\} =
		\mathbf\{v\}_r - \mathbf\{v\}_r\cdot\mathbf\{n\}\), \(\mathbf\{v\}_r =
		\mathbf\{v\}_j - \mathbf\{v\}_i\). The direction of the applied force
		is \(\mathbf\{t\} =
		\mathbf\{v\}_r - \mathbf\{v\}_r\cdot\mathbf\{n\}\{n\}\),
		\(\mathbf\{v\}_r = \mathbf\{v\}_j - \mathbf\{v\}_i\).
		The direction of the applied force is \(\mathbf\{t\} =
		\mathbf\{v_\{t,rel\}\}/\\|\mathbf\{v_\{t,rel\}\}\\|\) .

		The normal force value \(F_\{n0\}\) used to compute the critical force

src/GRANULAR/pair_granular.cpp

+8 −19

Original line number	Diff line number	Diff line
		@@ -305,6 +305,7 @@ void PairGranular::compute(int eflag, int vflag)

		delta = radsum - r;
		dR = delta*Reff;

		if (normal_model[itype][jtype] == JKR) {
		touch[jj] = 1;
		R2=Reff*Reff;
		@@ -1374,22 +1375,6 @@ double PairGranular::single(int i, int j, int itype, int jtype,
		vn2 = ny*vnnr;
		vn3 = nz*vnnr;

		double *rmass = atom->rmass;
		int *mask = atom->mask;
		mi = rmass[i];
		mj = rmass[j];
		if (fix_rigid) {
		if (mass_rigid[i] > 0.0) mi = mass_rigid[i];
		if (mass_rigid[j] > 0.0) mj = mass_rigid[j];
		}

		meff = mi*mj / (mi+mj);
		if (mask[i] & freeze_group_bit) meff = mj;
		if (mask[j] & freeze_group_bit) meff = mi;

		delta = radsum - r;
		dR = delta*Reff;

		// tangential component

		vt1 = vr1 - vn1;
		@@ -1407,6 +1392,9 @@ double PairGranular::single(int i, int j, int itype, int jtype,
		// if I or J part of rigid body, use body mass
		// if I or J is frozen, meff is other particle

		double *rmass = atom->rmass;
		int *mask = atom->mask;

		mi = rmass[i];
		mj = rmass[j];
		if (fix_rigid) {
		@@ -1557,7 +1545,8 @@ double PairGranular::single(int i, int j, int itype, int jtype,
		// rolling resistance
		//****************************************

		if (roll_model[itype][jtype] != ROLL_NONE) {
		if ((roll_model[itype][jtype] != ROLL_NONE)
		\|\| (twist_model[itype][jtype] != TWIST_NONE)) {
		relrot1 = omega[i][0] - omega[j][0];
		relrot2 = omega[i][1] - omega[j][1];
		relrot3 = omega[i][2] - omega[j][2];
		@@ -1623,7 +1612,7 @@ double PairGranular::single(int i, int j, int itype, int jtype,
		Mtcrit = mu_twist*Fncrit; // critical torque (eq 44)
		if (fabs(magtortwist) > Mtcrit) {
		magtortwist = -Mtcrit * signtwist; // eq 34
		}
		} else magtortwist = 0.0;
		}

		// set force and return no energy

Admin message