implement tangential force history in mindlin/force and mindlin_rescale/force (bf724332) · Commits · 郑智淋 / lammps

doc/src/pair_granular.rst

+96 −32

Original line number	Diff line number	Diff line
		@@ -140,7 +140,7 @@ where the force is computed as:

		\mathbf{F}_{ne, jkr} = \left(\frac{4Ea^3}{3R} - 2\pi a^2\sqrt{\frac{4\gamma E}{\pi a}}\right)\mathbf{n}

		Here, a is the radius of the contact zone, related to the overlap
		Here, :math:`a` is the radius of the contact zone, related to the overlap
		:math:`\delta` according to:

		.. math::
		@@ -167,7 +167,7 @@ following general form:

		\mathbf{F}_{n,damp} = -\eta_n \mathbf{v}_{n,rel}

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

		The optional damping keyword to the pair_coeff command followed by
		@@ -259,7 +259,9 @@ tangential model choices and their expected parameters are as follows:
		1. linear_nohistory : :math:`x_{\gamma,t}`, :math:`\mu_s`
		2. linear_history : :math:`k_t`, :math:`x_{\gamma,t}`, :math:`\mu_s`
		3. mindlin : :math:`k_t` or NULL, :math:`x_{\gamma,t}`, :math:`\mu_s`
		4. mindlin_rescale : :math:`k_t` or NULL, :math:`x_{\gamma,t}`, :math:`\mu_s`
		4. mindlin/force : :math:`k_t` or NULL, :math:`x_{\gamma,t}`, :math:`\mu_s`
		5. mindlin_rescale : :math:`k_t` or NULL, :math:`x_{\gamma,t}`, :math:`\mu_s`
		5. mindlin_rescale/force : :math:`k_t` or NULL, :math:`x_{\gamma,t}`, :math:`\mu_s`

		Here, :math:`x_{\gamma,t}` is a dimensionless multiplier for the normal
		damping :math:`\eta_n` that determines the magnitude of the tangential
		@@ -272,7 +274,7 @@ gran/hooke* style. The tangential force :math:`\mathbf{F}_t` is given by:

		.. math::

		\mathbf{F}_t = -min(\mu_t F_{n0}, \\|\mathbf{F}_\mathrm{t,damp}\\|) \mathbf{t}
		\mathbf{F}_t = -\min(\mu_t F_{n0}, \\|\mathbf{F}_\mathrm{t,damp}\\|) \mathbf{t}

		The tangential damping force :math:`\mathbf{F}_\mathrm{t,damp}` is given by:

		@@ -294,7 +296,7 @@ keyword also affects the tangential damping. The parameter
		literature use :math:`x_{\gamma,t} = 1` (:ref:`Marshall <Marshall2009>`,
		:ref:`Tsuji et al <Tsuji1992>`, :ref:`Silbert et al <Silbert2001>`). The relative
		tangential velocity at the point of contact is given by
		:math:`\mathbf{v}_{t, rel} = \mathbf{v}_{t} - (R_i\mathbf{\Omega}_i + R_j\mathbf{\Omega}_j) \times \mathbf{n}`, where :math:`\mathbf{v}_{t} = \mathbf{v}_r - \mathbf{v}_r\cdot\mathbf{n}\mathbf{n}`,
		:math:`\mathbf{v}_{t, rel} = \mathbf{v}_{t} - (R_i\mathbf{\Omega}_i + R_j\mathbf{\Omega}_j) \times \mathbf{n}`, where :math:`\mathbf{v}_{t} = \mathbf{v}_r - \mathbf{v}_r\cdot\mathbf{n}\ \mathbf{n}`,
		:math:`\mathbf{v}_r = \mathbf{v}_j - \mathbf{v}_i` .
		The direction of the applied force is :math:`\mathbf{t} = \mathbf{v_{t,rel}}/\\|\mathbf{v_{t,rel}}\\|` .

		@@ -314,22 +316,24 @@ form:

		.. math::

		F_{n0} = \\|\mathbf{F}_ne + 2 F_{pulloff}\\|
		F_{n0} = \\|\mathbf{F}_{ne} + 2 F_{pulloff}\\|

		Where :math:`F_{pulloff} = 3\pi \gamma R` for jkr\ , and
		:math:`F_{pulloff} = 4\pi \gamma R` for dmt\ .

		The remaining tangential options all use accumulated tangential displacement,
		or accumulated tangential force (i.e. contact history). This is discussed
		in details below for accumulated tangential displacement in the context
		of the linear_history option. The same treatment of the accumulated
		displacement, or accumulated force, applies to the other options as well.
		The remaining tangential options all use accumulated tangential
		displacement (i.e. contact history), except for the options
		mindlin/force and mindlin_rescale/force, that use accumulated
		tangential force instead, and are discussed further below.
		The accumulated tangential displacement is discussed in details below
		in the context of the linear_history option. The same treatment of
		the accumulated displacement applies to the other options as well.

		For tangential linear_history, the tangential force is given by:

		.. math::

		\mathbf{F}_t = -min(\mu_t F_{n0}, \\|-k_t\mathbf{\xi} + \mathbf{F}_\mathrm{t,damp}\\|) \mathbf{t}
		\mathbf{F}_t = -\min(\mu_t F_{n0}, \\|-k_t\mathbf{\xi} + \mathbf{F}_\mathrm{t,damp}\\|) \mathbf{t}

		Here, :math:`\mathbf{\xi}` is the tangential displacement accumulated
		during the entire duration of the contact:
		@@ -390,27 +394,18 @@ overlap region) to induce a torque on each particle according to:

		For tangential mindlin\ , the :ref:`Mindlin <Mindlin1949>` no-slip solution
		is used which differs from the linear_history option by an additional factor
		of a\ , the radius of the contact region. The tangential stiffness depends
		on the radius of the contact region and the force must therefore be computed
		incrementally. The accumulated tangential force characterizes the contact history.
		The increment of the elastic component of the tangential force :math:`\mathbf{F}_{te}`
		is given by:
		of :math:`a`, the radius of the contact region. The tangential force is given by:

		.. math::

		\mathrm{d}\mathbf{F}_{te} = -k_t a \mathbf{v}_{t,rel} \mathrm{d}\tau

		The tangential force is given by:

		.. math::
		\mathbf{F}_t = -\min(\mu_t F_{n0}, \\|-k_t a \mathbf{\xi} + \mathbf{F}_\mathrm{t,damp}\\|) \mathbf{t}

		\mathbf{F}_t = -min(\mu_t F_{n0}, \\|\mathbf{F}_{te} + \mathbf{F}_\mathrm{t,damp}\\|) \mathbf{t}

		Here, a is the radius of the contact region, given by :math:`a =\sqrt{R\delta}`
		Here, :math:`a` is the radius of the contact region, given by :math:`a =\sqrt{R\delta}`
		for all normal contact models, except for jkr\ , where it is given
		implicitly by :math:`\delta = a^2/R - 2\sqrt{\pi \gamma a/E}`, see
		discussion above. To match the Mindlin solution, one should set
		:math:`k_t = 8G_{eff}`, where :math:`G` is the effective shear modulus given by:
		:math:`k_t = 8G_{eff}`, where :math:`G_{eff}` is the effective shear modulus given by:

		.. math::

		@@ -424,23 +419,80 @@ normal contact model that specifies material parameters :math:`E` and
		case, mixing of the shear modulus for different particle types i and
		j is done according to the formula above.

		.. note::

		The radius of the contact region :math:`a` depends on the normal overlap.
		As a result, the tangential force for mindlin can change due to
		a variation in normal overlap, even with no change in tangential displacement.

		For tangential mindlin/force, the accumulated elastic tangential force
		characterizes the contact history, instead of the accumulated tangential
		displacement. This prevents the dependence of the tangential force on the
		normal overlap as noted above. The tangential force is given by:

		.. math::

		\mathbf{F}_t = -\min(\mu_t F_{n0}, \\|\mathbf{F}_{te} + \mathbf{F}_\mathrm{t,damp}\\|) \mathbf{t}

		The increment of the elastic component of the tangential force
		:math:`\mathbf{F}_{te}` is given by:

		.. math::

		\mathrm{d}\mathbf{F}_{te} = -k_t a \mathbf{v}_{t,rel} \mathrm{d}\tau

		The changes in frame of reference of the contacting pair of particles during
		contact are accounted for by the same formula as above, replacing the
		accumulated tangential displacement :math:`\xi`, by the accumulated tangential
		elastic force :math:`F_{te}`. When the tangential force exceeds the critical
		force, the tangential force is directly re-scaled to match the value for
		the critical force:

		.. math::

		\mathbf{F}_{te} = - \mu_t F_{n0}\mathbf{t} + \mathbf{F}_{t,damp}

		The same rules as those described for mindlin apply regarding the tangential
		stiffness and mixing of the shear modulus for different particle types.

		The mindlin_rescale option uses the same form as mindlin\ , but the
		magnitude of the tangential force is re-scaled as the contact
		magnitude of the tangential displacement is re-scaled as the contact
		unloads, i.e. if :math:`a < a_{t_{n-1}}`:

		.. math::

		\mathbf{F}_{te} = \mathbf{F}_{te, t_{n-1}} \frac{a}{a_{t_{n-1}}}
		\mathbf{\xi} = \mathbf{\xi_{t_{n-1}}} \frac{a}{a_{t_{n-1}}}

		Here, :math:`t_{n-1}` indicates the value at the previous time
		step. This rescaling accounts for the fact that a decrease in the
		contact area upon unloading leads to the contact being unable to
		support the previous tangential loading, and spurious energy is
		created without the rescaling above (:ref:`Walton <WaltonPC>` ). See also
		discussion in :ref:`Thornton et al, 2013 <Thornton2013>` , particularly
		equation 18(b) of that work and associated discussion.
		created without the rescaling above (:ref:`Walton <WaltonPC>` ).

		.. note::

		For mindlin, a decrease in the tangential force already occurs as the
		contact unloads, due to the dependence of the tangential force on the normal
		force described above. By re-scaling :math:`\xi`, mindlin_rescale
		effectively re-scales the tangential force twice, i.e., proportionally to
		:math:`a^2`. This peculiar behavior results from use of the accumulated
		tangential displacement to characterize the contact history. Although
		mindlin_rescale remains available for historic reasons and backward
		compatibility purposes, it should be avoided in favor of mindlin_rescale/force.

		The mindlin_rescale/force option uses the same form as mindlin/force,
		but the magnitude of the tangential elastic force is re-scaled as the contact
		unloads, i.e. if :math:`a < a_{t_{n-1}}`:

		.. math::

		\mathbf{F}_{te} = \mathbf{F}_{te, t_{n-1}} \frac{a}{a_{t_{n-1}}}

		This approach provides a better approximation of the :ref:`Mindlin-Deresiewicz <Mindlin1953>`
		laws and is more consistent than mindlin_rescale. See discussions in
		:ref:`Thornton et al, 2013 <Thornton2013>`, particularly equation 18(b) of that
		work and associated discussion, and :ref:`Agnolin and Roux, 2007 <AgnolinRoux2007>`,
		particularly Appendix A.

		----------

		@@ -477,7 +529,7 @@ exceeds a critical value:

		.. math::

		\mathbf{F}_{roll} = min(\mu_{roll} F_{n,0}, \\|\mathbf{F}_{roll,0}\\|)\mathbf{k}
		\mathbf{F}_{roll} = \min(\mu_{roll} F_{n,0}, \\|\mathbf{F}_{roll,0}\\|)\mathbf{k}

		Here, :math:`\mathbf{k} = \mathbf{v}_{roll}/\\|\mathbf{v}_{roll}\\|` is the direction of
		the pseudo-force. As with tangential displacement, the rolling
		@@ -529,7 +581,7 @@ is then truncated according to:

		.. math::

		\tau_{twist} = min(\mu_{twist} F_{n,0}, \tau_{twist,0})
		\tau_{twist} = \min(\mu_{twist} F_{n,0}, \tau_{twist,0})

		Similar to the sliding and rolling displacement, the angular
		displacement is rescaled so that it corresponds to the critical value
		@@ -794,3 +846,15 @@ Technology, 233, 30-46.
		.. _WaltonPC:

		(Otis R. Walton) Walton, O.R., Personal Communication

		** _Mindlin1953:

		(Mindlin and Deresiewicz, 1953) Mindlin, R.D., & Deresiewicz, H (1953).
		Elastic Spheres in Contact under Varying Oblique Force.
		J. Appl. Mech., ASME 20, 327-344.

		** _AgnolinRoux2007:

		(Agnolin and Roux 2007) Agnolin, I. & Roux, J-N. (2007).
		Internal states of model isotropic granular packings.
		I. Assembling process, geometry, and contact networks. Phys. Rev. E, 76, 061302.
		No newline at end of file

src/GRANULAR/pair_granular.cpp

+59 −24

Original line number	Diff line number	Diff line
		@@ -54,7 +54,8 @@ using namespace MathSpecial;
		enum {HOOKE, HERTZ, HERTZ_MATERIAL, DMT, JKR};
		enum {VELOCITY, MASS_VELOCITY, VISCOELASTIC, TSUJI};
		enum {TANGENTIAL_NOHISTORY, TANGENTIAL_HISTORY,
		TANGENTIAL_MINDLIN, TANGENTIAL_MINDLIN_RESCALE};
		TANGENTIAL_MINDLIN, TANGENTIAL_MINDLIN_RESCALE,
		TANGENTIAL_MINDLIN_FORCE, TANGENTIAL_MINDLIN_RESCALE_FORCE};
		enum {TWIST_NONE, TWIST_SDS, TWIST_MARSHALL};
		enum {ROLL_NONE, ROLL_SDS};

		@@ -377,8 +378,11 @@ void PairGranular::compute(int eflag, int vflag)
		// tangential force, including history effects
		//****************************************

		// For linear forces, history = cumulative tangential displacement
		// For nonlinear forces, history = cumulative elastic tangential force
		// For linear, mindlin, mindlin_rescale:
		// history = cumulative tangential displacement
		//
		// For mindlin/force, mindlin_rescale/force:
		// history = cumulative tangential elastic force

		// tangential component
		vt1 = vr1 - vn1;
		@@ -422,12 +426,15 @@ void PairGranular::compute(int eflag, int vflag)
		damp_normal_prefactor;

		if (tangential_history) {
		if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN) {
		if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN \|\|
		tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_FORCE) {
		k_tangential *= a;
		} else if (tangential_model[itype][jtype] ==
		TANGENTIAL_MINDLIN_RESCALE) {
		TANGENTIAL_MINDLIN_RESCALE \|\|
		tangential_model[itype][jtype] ==
		TANGENTIAL_MINDLIN_RESCALE_FORCE) {
		k_tangential *= a;
		// on unloading, rescale the shear force
		// on unloading, rescale the shear displacements/force
		if (a < history[3]) {
		double factor = a/history[3];
		history[0] *= factor;
		@@ -435,11 +442,11 @@ void PairGranular::compute(int eflag, int vflag)
		history[2] *= factor;
		}
		}
		// rotate and update displacements / force
		// rotate and update displacements / force.
		// see e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235
		if (historyupdate) {
		rsht = history[0]nx + history[1]ny + history[2]*nz;
		if (fabs(rsht) < EPSILON) rsht = 0; // EPSILON still valid when history is a force???
		if (fabs(rsht) < EPSILON) rsht = 0; // EPSILON = 1e-10 can fail for small granular media:
		if (rsht > 0) {
		shrmag = sqrt(history[0]history[0] + history[1]history[1] +
		history[2]*history[2]);
		@@ -454,22 +461,31 @@ void PairGranular::compute(int eflag, int vflag)
		history[1] *= scalefac;
		history[2] *= scalefac;
		}
		// update tangential displacement / force history
		if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY) {
		// update history
		if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY \|\|
		tangential_model[itype][jtype] == TANGENTIAL_MINDLIN \|\|
		tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) {
		// tangential displacement
		history[0] += vtr1*dt;
		history[1] += vtr2*dt;
		history[2] += vtr3*dt;
		} else { // Eq. (18) Thornton et al., 2013 (http://dx.doi.org/10.1016/j.powtec.2012.08.012)
		} else {
		// tangential force
		// see e.g. eq. 18 of Thornton et al, Pow. Tech. 2013, v223,p30-46
		history[0] -= k_tangentialvtr1dt;
		history[1] -= k_tangentialvtr2dt;
		history[2] -= k_tangentialvtr3dt;
		if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE)
		history[3] = a;
		}
		if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE \|\|
		tangential_model[itype][jtype] ==
		TANGENTIAL_MINDLIN_RESCALE_FORCE)
		history[3] = a;
		}

		// tangential forces = history + tangential velocity damping
		if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY) {
		if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY \|\|
		tangential_model[itype][jtype] == TANGENTIAL_MINDLIN \|\|
		tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) {
		fs1 = -k_tangentialhistory[0] - damp_tangentialvtr1;
		fs2 = -k_tangentialhistory[1] - damp_tangentialvtr2;
		fs3 = -k_tangentialhistory[2] - damp_tangentialvtr3;
		@@ -485,7 +501,10 @@ void PairGranular::compute(int eflag, int vflag)
		shrmag = sqrt(history[0]history[0] + history[1]history[1] +
		history[2]*history[2]);
		if (shrmag != 0.0) {
		if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY) {
		if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY \|\|
		tangential_model[itype][jtype] == TANGENTIAL_MINDLIN \|\|
		tangential_model[itype][jtype] ==
		TANGENTIAL_MINDLIN_RESCALE) {
		history[0] = -1.0/k_tangential(Fscritfs1/fs +
		damp_tangential*vtr1);
		history[1] = -1.0/k_tangential(Fscritfs2/fs +
		@@ -760,7 +779,8 @@ void PairGranular::coeff(int narg, char **arg)

		//Defaults
		normal_model_one = tangential_model_one = -1;
		roll_model_one = twist_model_one = 0;
		roll_model_one = ROLL_NONE;
		twist_model_one = TWIST_NONE;
		damping_model_one = VISCOELASTIC;

		int iarg = 2;
		@@ -846,7 +866,9 @@ void PairGranular::coeff(int narg, char **arg)
		iarg += 4;
		} else if ((strcmp(arg[iarg+1], "linear_history") == 0) \|\|
		(strcmp(arg[iarg+1], "mindlin") == 0) \|\|
		(strcmp(arg[iarg+1], "mindlin_rescale") == 0)) {
		(strcmp(arg[iarg+1], "mindlin_rescale") == 0) \|\|
		(strcmp(arg[iarg+1], "mindlin/force") == 0) \|\|
		(strcmp(arg[iarg+1], "mindlin_rescale/force") == 0)) {
		if (iarg + 4 >= narg)
		error->all(FLERR,"Illegal pair_coeff command, "
		"not enough parameters provided for tangential model");
		@@ -856,9 +878,15 @@ void PairGranular::coeff(int narg, char **arg)
		tangential_model_one = TANGENTIAL_MINDLIN;
		else if (strcmp(arg[iarg+1], "mindlin_rescale") == 0)
		tangential_model_one = TANGENTIAL_MINDLIN_RESCALE;
		else if (strcmp(arg[iarg+1], "mindlin/force") == 0)
		tangential_model_one = TANGENTIAL_MINDLIN_FORCE;
		else if (strcmp(arg[iarg+1], "mindlin_rescale/force") == 0)
		tangential_model_one = TANGENTIAL_MINDLIN_RESCALE_FORCE;
		tangential_history = 1;
		if ((tangential_model_one == TANGENTIAL_MINDLIN \|\|
		tangential_model_one == TANGENTIAL_MINDLIN_RESCALE) &&
		tangential_model_one == TANGENTIAL_MINDLIN_RESCALE \|\|
		tangential_model_one == TANGENTIAL_MINDLIN_FORCE \|\|
		tangential_model_one == TANGENTIAL_MINDLIN_RESCALE_FORCE) &&
		(strcmp(arg[iarg+2], "NULL") == 0)) {
		if (normal_model_one == HERTZ \|\| normal_model_one == HOOKE) {
		error->all(FLERR, "NULL setting for Mindlin tangential "
		@@ -1040,7 +1068,8 @@ void PairGranular::init_style()
		}
		for (int i = 1; i <= atom->ntypes; i++)
		for (int j = i; j <= atom->ntypes; j++)
		if (tangential_model[i][j] == TANGENTIAL_MINDLIN_RESCALE) {
		if (tangential_model[i][j] == TANGENTIAL_MINDLIN_RESCALE \|\|
		tangential_model[i][j] == TANGENTIAL_MINDLIN_RESCALE_FORCE) {
		size_history += 1;
		roll_history_index += 1;
		twist_history_index += 1;
		@@ -1510,6 +1539,12 @@ double PairGranular::single(int i, int j, int itype, int jtype,
		// tangential force, including history effects
		//****************************************

		// For linear, mindlin, mindlin_rescale:
		// history = cumulative tangential displacement
		//
		// For mindlin/force, mindlin_rescale/force:
		// history = cumulative tangential elastic force

		// tangential component
		vt1 = vr1 - vn1;
		vt2 = vr2 - vn2;
		@@ -1545,9 +1580,7 @@ double PairGranular::single(int i, int j, int itype, int jtype,
		damp_tangential = tangential_coeffs[itype][jtype][1]*damp_normal_prefactor;

		if (tangential_history) {
		if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN) {
		k_tangential *= a;
		} else if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) {
		if (tangential_model[itype][jtype] != TANGENTIAL_HISTORY) {
		k_tangential *= a;
		}

		@@ -1555,7 +1588,9 @@ double PairGranular::single(int i, int j, int itype, int jtype,
		history[2]*history[2]);

		// tangential forces = history + tangential velocity damping
		if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY) {
		if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY \|\|
		tangential_model[itype][jtype] == TANGENTIAL_MINDLIN \|\|
		tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) {
		fs1 = -k_tangentialhistory[0] - damp_tangentialvtr1;
		fs2 = -k_tangentialhistory[1] - damp_tangentialvtr2;
		fs3 = -k_tangentialhistory[2] - damp_tangentialvtr3;

Admin message