Commit cc14103f authored by Jibril B. Coulibaly's avatar Jibril B. Coulibaly
Browse files

Bug fixes in granular pair style: 

- correct formula for tangent forces in style with no history in compute() and in single() functions
- remove tangent history update in the single() function
- implement correct output for tangent, normal and rolling forces in single() function
- correct typos in documentation
parent 37a046cf

doc/src/pair_granular.txt

+4 −4
Original line number Diff line number Diff line
@@ -100,7 +100,7 @@ on particle {i} due to contact with particle {j} is given by: 
\mathbf\{F\}_\{ne, Hooke\} = k_N \delta_\{ij\} \mathbf\{n\}

\end\{equation\}



Where \(\delta = R_i + R_j - \|\mathbf\{r\}_\{ij\}\|\) is the particle

Where \(\delta_\{ij\} = R_i + R_j - \|\mathbf\{r\}_\{ij\}\|\) is the particle

overlap, \(R_i, R_j\) are the particle radii, \(\mathbf\{r\}_\{ij\} =

\mathbf\{r\}_i - \mathbf\{r\}_j\) is the vector separating the two

particle centers (note the i-j ordering so that \(F_\{ne\}\) is
@@ -411,8 +411,8 @@ option by an additional factor of {a}, the radius of the contact region. The tan 
\mathbf\{F\}_t =  -min(\mu_t F_\{n0\}, \|-k_t a \mathbf\{\xi\} + \mathbf\{F\}_\mathrm\{t,damp\}\|) \mathbf\{t\}

\end\{equation\}



Here, {a} is the radius of the contact region, given by \(a = \delta

R\) for all normal contact models, except for {jkr}, where it is given

Here, {a} is the radius of the contact region, given by \(a =\sqrt\{R\delta\}\)

 for all normal contact models, except for {jkr}, where it is given

implicitly by \(\delta = a^2/R - 2\sqrt\{\pi \gamma a/E\}\), see

discussion above. To match the Mindlin solution, one should set \(k_t

= 8G\), where \(G\) is the shear modulus, related to Young's modulus
@@ -680,7 +680,7 @@ The single() function of these pair styles returns 0.0 for the energy 
of a pairwise interaction, since energy is not conserved in these

dissipative potentials.  It also returns only the normal component of

the pairwise interaction force.  However, the single() function also

calculates 10 extra pairwise quantities.  The first 3 are the

calculates 12 extra pairwise quantities.  The first 3 are the

components of the tangential force between particles I and J, acting

on particle I.  The 4th is the magnitude of this tangential force.

The next 3 (5-7) are the components of the rolling torque acting on

src/GRANULAR/pair_granular.cpp

+19 −23
Original line number Diff line number Diff line
@@ -391,6 +391,7 @@ void PairGranular::compute(int eflag, int vflag) 
        } else {

          Fncrit = fabs(Fntot);

        }

		Fscrit = tangential_coeffs[itype][jtype][2] * Fncrit;



        //------------------------------

        // tangential forces
@@ -446,7 +447,6 @@ void PairGranular::compute(int eflag, int vflag) 
          fs3 = -k_tangential*history[2] - damp_tangential*vtr3;



          // rescale frictional displacements and forces if needed

          Fscrit = tangential_coeffs[itype][jtype][2] * Fncrit;

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

          if (fs > Fscrit) {

            shrmag = sqrt(history[0]*history[0] + history[1]*history[1] +
@@ -464,8 +464,8 @@ void PairGranular::compute(int eflag, int vflag) 
            } else fs1 = fs2 = fs3 = 0.0;

          }

        } else { // classic pair gran/hooke (no history)

          fs = meff*damp_tangential*vrel;

          if (vrel != 0.0) Ft = MIN(Fne,fs) / vrel;

          fs = damp_tangential*vrel; // From documentation: F_{t,damp} = - \eta_t v_{t,rel}, no need for extra `meff`

          if (vrel != 0.0) Ft = MIN(Fscrit,fs) / vrel; // From documentation: critical force `Fscrit` used, not elastic normal force `Fne`

          else Ft = 0.0;

          fs1 = -Ft*vtr1;

          fs2 = -Ft*vtr2;
@@ -635,7 +635,7 @@ void PairGranular::compute(int eflag, int vflag) 
            torque[j][2] -= torroll3;

          }

        }

        if (evflag) ev_tally_xyz(i,j,nlocal,0,

        if (evflag) ev_tally_xyz(i,j,nlocal,0,//Should `newton_pair` passed instead of 0 ?

            0.0,0.0,fx,fy,fz,delx,dely,delz);

      }

    }
@@ -1451,11 +1451,13 @@ double PairGranular::single(int i, int j, int itype, int jtype, 
  }



  if (damping_model[itype][jtype] == VELOCITY) {

    damp_normal = normal_coeffs[itype][jtype][1];

    damp_normal = 1;

  } else if (damping_model[itype][jtype] == MASS_VELOCITY) {

    damp_normal = meff;

  } else if (damping_model[itype][jtype] == VISCOELASTIC) {

    damp_normal = normal_coeffs[itype][jtype][1]*a*meff;

    damp_normal = a*meff;

  } else if (damping_model[itype][jtype] == TSUJI) {

    damp_normal = normal_coeffs[itype][jtype][1]*sqrt(meff*knfac);

    damp_normal = sqrt(meff*knfac);

  }



  damp_normal_prefactor = normal_coeffs[itype][jtype][1]*damp_normal;
@@ -1473,6 +1475,7 @@ double PairGranular::single(int i, int j, int itype, int jtype, 
      if (neighprev >= jnum) neighprev = 0;

      if (jlist[neighprev] == j) break;

    }

	// the `history` pointer must not be modified here in single() function. already calculated in the compute() function. If modified here it changes the pair forces that have friction/twisting/rolling and history effects !

    history = &allhistory[size_history*neighprev];

  }
@@ -1506,6 +1509,7 @@ double PairGranular::single(int i, int j, int itype, int jtype, 
  } else {

    Fncrit = fabs(Fntot);

  }

  Fscrit = tangential_coeffs[itype][jtype][2] * Fncrit;



  //------------------------------

  // tangential forces
@@ -1518,13 +1522,6 @@ double PairGranular::single(int i, int j, int itype, int jtype, 
      k_tangential *= a;

    } else if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) {

      k_tangential *= a;

      // on unloading, rescale the shear displacements

      if (a < history[3]) {

        double factor = a/history[3];

        history[0] *= factor;

        history[1] *= factor;

        history[2] *= factor;

      }

    }



    shrmag = sqrt(history[0]*history[0] + history[1]*history[1] +
@@ -1535,28 +1532,26 @@ double PairGranular::single(int i, int j, int itype, int jtype, 
    fs2 = -k_tangential*history[1] - damp_tangential*vtr2;

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



    // rescale frictional displacements and forces if needed

    Fscrit = tangential_coeffs[itype][jtype][2] * Fncrit;

    // rescale frictional forces if needed

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

    if (fs > Fscrit) {

      if (shrmag != 0.0) {

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

        fs1 *= Fscrit/fs;

        fs2 *= Fscrit/fs;

        fs3 *= Fscrit/fs;

      } else fs1 = fs2 = fs3 = 0.0;

		fs *= Fscrit/fs; // saves the correct value of `fs` to svector

      } else fs1 = fs2 = fs3 = fs = 0.0; // saves the correct of `fs` value to svector

    }



  // classic pair gran/hooke (no history)

  } else {

    fs = meff*damp_tangential*vrel;

    if (vrel != 0.0) Ft = MIN(Fne,fs) / vrel;

    fs = damp_tangential*vrel;

    if (vrel != 0.0) Ft = MIN(Fscrit,fs) / vrel;

    else Ft = 0.0;

    fs1 = -Ft*vtr1;

    fs2 = -Ft*vtr2;

    fs3 = -Ft*vtr3;

	fs = Ft*vrel; // saves the correct value of `fs` to svector

  }



  //****************************************
@@ -1601,7 +1596,8 @@ double PairGranular::single(int i, int j, int itype, int jtype, 
        fr1 *= Frcrit/fr;

        fr2 *= Frcrit/fr;

        fr3 *= Frcrit/fr;

      } else fr1 = fr2 = fr3 = 0.0;

		fr *= Frcrit/fr; // saves the correct value of `fr` to svector

      } else fr1 = fr2 = fr3 = fr = 0.0; // saves the correct value of `fr` to svector

    }



  }

CRA Git | Maintained and supported by SUSTech CRA and CCSE