Rolling friction torque

diego.peinado's picture
Submitted by diego.peinado on Mon, 01/09/2012 - 18:04

Hello, I found some strange behavior regarding the way the rolling friction torque is calculated.
I use a lattice command to create the particles. Then they are in contact, but there is no shear at all in the first step.
This is the code used to create the particles:

lattice fcc 0.01414213562
region cil cylinder y 0 0.01 0.79 0 0.1 units box
region res prism -1.1 1.1 0.01 0.79 -1.1 -0.20 0.0 0.0 0.0 units box
region part intersect 2 cil res units box
create_atoms 1 region part units box
group myPart region part
set group myPart diameter 0.01
set group myPart density 2500

The current code is:
if(rollingflag)
{
reff=radi*radj/(radi+radj);
wrmag = sqrt(wr1*wr1+wr2*wr2+wr3*wr3);
if (wrmag > 0.)
{
r_torque[0] = rmu*kn*deltan*wr1/wrmag*reff;
r_torque[1] = rmu*kn*deltan*wr2/wrmag*reff;
r_torque[2] = rmu*kn*deltan*wr3/wrmag*reff;
// remove normal (torsion) part of torque
r_torque_n[0] = r_torque[0] * delx * rinv;
r_torque_n[1] = r_torque[1] * dely * rinv;
r_torque_n[2] = r_torque[2] * delz * rinv;
vectorSubtract3D(r_torque,r_torque_n,r_torque);
}
}
else vectorZeroize3D(r_torque);
In this case, rollingflat is true, but wrmag = 0. So there is no inicialization to r_torque[]. Then, when the following code is executed
torque[i][0] -= cri*tor1 + r_torque[0];
torque[i][1] -= cri*tor2 + r_torque[1];
torque[i][2] -= cri*tor3 + r_torque[2];
There can be any value in r_torque, giving numerical errors in this evaluation (inf, nan, even segmentation fault).
If the code is changed as:

vectorZeroize3D(r_torque);
if(rollingflag)
{
reff=radi*radj/(radi+radj);
wrmag = sqrt(wr1*wr1+wr2*wr2+wr3*wr3);
if (wrmag > 0.)
{
r_torque[0] = rmu*kn*deltan*wr1/wrmag*reff;
r_torque[1] = rmu*kn*deltan*wr2/wrmag*reff;
r_torque[2] = rmu*kn*deltan*wr3/wrmag*reff;
// remove normal (torsion) part of torque
r_torque_n[0] = r_torque[0] * delx * rinv;
r_torque_n[1] = r_torque[1] * dely * rinv;
r_torque_n[2] = r_torque[2] * delz * rinv;
vectorSubtract3D(r_torque,r_torque_n,r_torque);
}
}
Then always is inicializated. I've only found this place to correct, but perhaps are more.

Thanks for your efforts, cheers

ckloss's picture

ckloss | Wed, 01/11/2012 - 22:24

Hello Diego,

hope you had happy holidays and a good start into the new year!
Thanks for the info - yes, you are correct, I will fix that!

Cheers, Christoph