Discussion: Implementation of tangential spring in pair_style gran/hooke

Submitted by Yansan on Mon, 11/08/2010 - 17:28

Hi all,
I think there are some problems with the implementation of the tangential spring in the pair_style gran/hooke, especially for soft particles where the overlap between the particles gets high.
So far the relative rotational velocities between 2 particles are calculated by using the radius of the particles:


wr1 = (radi*omega[i][0] + radj*omega[j][0]) * rinv;

The force application point for the tangential forces also lies on the radius of the particles:

torque[i][0] -= radi*tor1;

In my opinion the relative rotational velocities should be calculated by using the distance between the contact point of the two particles and the centre point of each particle. The Formula should look like:

wr1 = ((radi-0.5*delta)*omega[i][0] + (radj-0.5*delta)*omega[j][0]) * rinv;

Where delta is the overlap between particle i and j (assuming that the contact point is in the middle of the overlap). This applies to the force application point of the tangential forces too.
In my opinion the tangential forces should apply at the contact point and not at the radius of the particles.
What do you think about this?

King regards
Yansan

Forums:

ckloss's picture

ckloss | Mon, 11/08/2010 - 17:38

>>for soft particles where the overlap between the particles gets high.
Normally it is recommended to have an inter penetration of ~0.5% - so the error should be small. But I agree - yI would also say that this would probably be better. There is also some "incorrectness" in the zero-order rotation of the shear history - but it is computationally cheap.

If you are interested in this effect, it would be cool if you would try to compare these two implementations so we have some facts to discuss about.

Cheers & all the best,
Christoph

marketos | Mon, 11/15/2010 - 16:45

Hi Yansan,

I agree with you that this is an issue with the implementation. As Christoph said the difference should be clearly visible for large interpenetrations and have a negligible effect for small interpenetrations. I think that the code definitely needs to be modified, as in our simulations we have cases where the forces on the grains (and so the overlap) will be large. I am currently looking at how to do this- it should be roughly as you said, subtracting overlap/2 from each radius and using this for calculating the relative tangential velocities due to rotations and the torques created by tangential forces.

We need to be careful so as to minimise the number of computations while doing this - this part of the calculation is run for each particle pair every timestep so inefficiencies here will be significant but unfortunately I am not an expert on this. I will let you know when I have something that works and please do the same.

Thanks,

George

Yansan | Mon, 11/15/2010 - 18:02

Hi George,
I am not an expert too but I do it like this:


//declaration at the beginning of the compute
double cri, crj;
double deltan; //originally in line 252 (in pair_gran_hooke_history.cpp)
[...]
//if particles are in contact starting from line 202
deltan=radsum-r;
[...]
// relative rotational velocity
cri = radi-0.5*deltan;
crj = radj-0.5*deltan;
wr1 = (cri*omega[i][0] + crj*omega[j][0]) * rinv;
wr2 = (cri*omega[i][1] + crj*omega[j][1]) * rinv;
wr3 = (cri*omega[i][2] + crj*omega[j][2]) * rinv;
[...]
tor1 = rinv * (dely*fs3 - delz*fs2);
tor2 = rinv * (delz*fs1 - delx*fs3);
tor3 = rinv * (delx*fs2 - dely*fs1);
torque[i][0] -= cri*tor1;
torque[i][1] -= cri*tor2;
torque[i][2] -= cri*tor3;

if (j < nlocal) {
f[j][0] -= fx;
f[j][1] -= fy;
f[j][2] -= fz;
torque[j][0] -= crj*tor1;
torque[j][1] -= crj*tor2;
torque[j][2] -= crj*tor3;
}

This works fine for me and in summary you have 2 more variables and a 2 cheap calculations. I hope you are able to follow my instructions.

Kind regards
Yansan

ckloss's picture

ckloss | Tue, 11/16/2010 - 10:07

Hi Yansan,

thanks for sharing your findings - I would like your version and the version in the main branch. We can do that as soon as the testharness is up and running on a bigger scale.

Cheers,
Christoph

marketos | Tue, 11/16/2010 - 16:39

Hi Yansan,

Thank you for the code you sent. I think this is the way forward but I agree with Christoph, we need to be 100% sure that the equations are fine before he makes a change to the version he releases. (But I think they are fine)

I assume that you are using walls for your simulations. In this case there is another compute function inside the file fix_wall_gran_hooke_history.cpp . This needs to be modified too, as again here radius instead of the distance between the wall and the ball (r) is used for the relevant calculations. I spent most of today checking what the changes should be and I am 99.9% sure that you need to change the following...

In file fix_wall_gran_hooke_history.cpp:
Comment out these lines- they are not needed any more as r should be used instead of radius, making wr1=omega[0] etc..
//wr1 = radius*omega[0] * rinv;
//wr2 = radius*omega[1] * rinv;
//wr3 = radius*omega[2] * rinv;

Then change the lines where vtr1, vtr2, vtr3 are calculated
vtr1 = vt1 - dz*omega[1]+dy*omega[2];
vtr2 = vt2 - dx*omega[2]+dz*omega[0];
vtr3 = vt3 - dy*omega[0]+dx*omega[1];

Do something similar for the torques- comment out the lines below as r instead of radius should be used when updating torques and rinv*r=1
//tor1 = rinv * (dy*fs3 - dz*fs2);
//tor2 = rinv * (dz*fs1 - dx*fs3);
//tor3 = rinv * (dx*fs2 - dy*fs1);

Then modify the lines where torque is calculated:

torque[0] -= dy*fs3 - dz*fs2;
torque[1] -= dz*fs1 - dx*fs3;
torque[2] -= dx*fs2 - dy*fs1;

And finally delete the definition of variables not needed any more

double meff,damp,ccel,vtr1,vtr2,vtr3,vrel;//wr1,wr2,wr3,meff,damp,ccel,vtr1,vtr2,vtr3,vrel;
double fn,fs,fs1,fs2,fs3,fx,fy,fz;//,tor1,tor2,tor3;

I would be very grateful if you could check to see if I have made any mistake. I think the formulae are OK.

One final thing to note is that I do not know what happens when the wall rotates- presumably this is taken into consideration somewhere else and modifications have to be made there, or the angular velocity of the wall needs to be included in these formulae. It is also possible that other files need to be modified accordingly- Christoph any ideas?

Thanks,

George

ckloss's picture

ckloss | Tue, 11/16/2010 - 16:49

>>It is also possible that other files need to be modified accordingly- Christoph any ideas?
I don't think so!

Thanks you all for this valuable discussion!

Cheers,
Christoph

marketos | Tue, 11/16/2010 - 18:19

I guess files pair_gran_hooke.cpp and fix_wall_gran_hooke.cpp would also need to be modified for consistency, but I am not sure if anybody uses them at the moment. I cannot think of anything else at the moment.

Thanks for all your help,

George

Yansan | Wed, 11/17/2010 - 09:44

Hi George,
Of course you are right and you have to change fix_wall_gran_hooke_history.cpp as well, I just forgot to mention it. I changed my code in a similar way but your code looks fine for me and should be faster then mine (I did it like in pair_gran_hooke_history.cpp).

Thanks,
Yansan