Replace accumulated displacement by accumulated force for tangential force in...

.. math::

   \mathbf{F}_{ne, Hooke} = k_N \delta_{ij} \mathbf{n}

   \mathbf{F}_{ne, Hooke} = k_n \delta_{ij} \mathbf{n}

Where :math:`\delta_{ij} = R_i + R_j - \|\mathbf{r}_{ij}\|` is the particle

overlap, :math:`R_i, R_j` are the particle radii, :math:`\mathbf{r}_{ij} = \mathbf{r}_i - \mathbf{r}_j` is the vector separating the two

.. math::

   \mathbf{F}_{ne, Hertz} = k_N R_{eff}^{1/2}\delta_{ij}^{3/2} \mathbf{n}

   \mathbf{F}_{ne, Hertz} = k_n R_{eff}^{1/2}\delta_{ij}^{3/2} \mathbf{n}

Here, :math:`R_{eff} = \frac{R_i R_j}{R_i + R_j}` is the effective

radius, denoted for simplicity as *R* from here on.  For *hertz*\ , the

modulus, with :math:`\nu_i, \nu_j` the Poisson ratios of the particles of

types *i* and *j*\ . Note that if the elastic modulus and the shear

modulus of the two particles are the same, the *hertz/material* model

is equivalent to the *hertz* model with :math:`k_N = 4/3 E_{eff}`

is equivalent to the *hertz* model with :math:`k_n = 4/3 E_{eff}`

The *dmt* model corresponds to the

:ref:`(Derjaguin-Muller-Toporov) <DMT1975>` cohesive model, where the force

For *tangential linear_nohistory*, a simple velocity-dependent Coulomb

friction criterion is used, which mimics the behavior of the *pair

gran/hooke* style. The tangential force (\mathbf{F}_t\) is given by:

gran/hooke* style. The tangential force :math:`\mathbf{F}_t` is given by:

.. math::

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\Omega_i + R_j\Omega_j) \times \mathbf{n}`, where :math:`\mathbf{v}_{t} = \mathbf{v}_r - \mathbf{v}_r\cdot\mathbf{n}{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}}\|` .

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 (i.e. contact history). This is discussed below in the

context of the *linear_history* option, but the same treatment of the

accumulated displacement applies to the other options as well.

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.

For *tangential linear_history*, the tangential force is given by:

.. math::

   \mathbf{\xi} = -\frac{1}{k_t}\left(\mu_t F_{n0}\mathbf{t} + \mathbf{F}_{t,damp}\right)

   \mathbf{\xi} = -\frac{1}{k_t}\left(\mu_t F_{n0}\mathbf{t} - \mathbf{F}_{t,damp}\right)

The tangential force is added to the total normal force (elastic plus

damping) to produce the total force on the particle. The tangential

   \mathbf{\tau}_j = -(R_j - 0.5 \delta) \mathbf{n} \times \mathbf{F}_t

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 force is given by:

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:

.. math::

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

   \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}, \|\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}`

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 = 4G/(2-\nu)`, where :math:`G` is the shear modulus, related to Young's modulus

:math:`E` by :math:`G = E/(2(1+\nu))`, where :math:`\nu` is Poisson's ratio. This

can also be achieved by specifying *NULL* for :math:`k_t`, in which case a

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::

   G_{eff} = \left(\frac{2-\nu_i}{G_i} + \frac{2-\nu_j}{G_j}\right)^{-1}

where :math:`G` is the shear modulus, related to Young's modulus :math:`E`

and Poisson's ratio :math:`\nu` by :math:`G = E/(2(1+\nu))`. This can also be

achieved by specifying *NULL* for :math:`k_t`, in which case a

normal contact model that specifies material parameters :math:`E` and

:math:`\nu` is required (e.g. *hertz/material*\ , *dmt* or *jkr*\ ). In this

case, mixing of the shear modulus for different particle types *i* and

*j* is done according to:

*j* is done according to the formula above.

.. math::

   1/G = 2(2-\nu_i)(1+\nu_i)/E_i + 2(2-\nu_j)(1+\nu_j)/E_j

The *mindlin_rescale* option uses the same form as *mindlin*\ , but the

magnitude of the tangential displacement is re-scaled as the contact

magnitude of the tangential force is re-scaled as the contact

unloads, i.e. if :math:`a < a_{t_{n-1}}`:

.. math::

   \mathbf{\xi} = \mathbf{\xi_{t_{n-1}}} \frac{a}{a_{t_{n-1}}}

   \mathbf{F}_{te} = \mathbf{F}_{te, 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

        // tangential force, including history effects

        //****************************************

        // For linear forces, history = cumulative tangential displacement

        // For nonlinear forces, history = cumulative elastic tangential force

        // tangential component

        vt1 = vr1 - vn1;

        vt2 = vr2 - vn2;

          } else if (tangential_model[itype][jtype] ==

                     TANGENTIAL_MINDLIN_RESCALE) {

            k_tangential *= a;

            // on unloading, rescale the shear displacements

            // on unloading, rescale the shear force

            if (a < history[3]) {

              double factor = a/history[3];

              history[0] *= factor;

              history[2] *= factor;

            }

          }

          // rotate and update displacements.

          // 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;

            if (fabs(rsht) < EPSILON) rsht = 0; // EPSILON still valid when history is a force???

            if (rsht > 0) {

              shrmag = sqrt(history[0]*history[0] + history[1]*history[1] +

                                               history[2]*history[2]);

              history[1] *= scalefac;

              history[2] *= scalefac;

            }

            // update history

            // update tangential displacement / force history

            if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY) {

              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)

              history[0] -= k_tangential*vtr1*dt;

              history[1] -= k_tangential*vtr2*dt;

              history[2] -= k_tangential*vtr3*dt;

              if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE)

                history[3] = a;

            }

          }

          // tangential forces = history + tangential velocity damping

          if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY) {

            fs1 = -k_tangential*history[0] - damp_tangential*vtr1;

            fs2 = -k_tangential*history[1] - damp_tangential*vtr2;

            fs3 = -k_tangential*history[2] - damp_tangential*vtr3;

          } else {

            fs1 = history[0] - damp_tangential*vtr1;

            fs2 = history[1] - damp_tangential*vtr2;

            fs3 = history[2] - damp_tangential*vtr3;

          }

          // rescale frictional displacements and forces if needed

          fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3);

            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) {

                history[0] = -1.0/k_tangential*(Fscrit*fs1/fs +

                                                damp_tangential*vtr1);

                history[1] = -1.0/k_tangential*(Fscrit*fs2/fs +

                                                damp_tangential*vtr2);

                history[2] = -1.0/k_tangential*(Fscrit*fs3/fs +

                                                damp_tangential*vtr3);

              } else {

                history[0] = Fscrit*fs1/fs + damp_tangential*vtr1;

                history[1] = Fscrit*fs2/fs + damp_tangential*vtr2;

                history[2] = Fscrit*fs3/fs + damp_tangential*vtr3;

              }

              fs1 *= Fscrit/fs;

              fs2 *= Fscrit/fs;

              fs3 *= Fscrit/fs;

        history[2]*history[2]);

    // tangential forces = history + tangential velocity damping

    if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY) {

      fs1 = -k_tangential*history[0] - damp_tangential*vtr1;

      fs2 = -k_tangential*history[1] - damp_tangential*vtr2;

      fs3 = -k_tangential*history[2] - damp_tangential*vtr3;

    } else {

      fs1 = history[0] - damp_tangential*vtr1;

      fs2 = history[1] - damp_tangential*vtr2;

      fs3 = history[2] - damp_tangential*vtr3;

    }

    // rescale frictional forces if needed

    fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3);