← Back to team overview

yade-users team mailing list archive

Re: Contact Model and Contact Stiffness


mohammad islam a écrit :
 Dear All,
I am losing myself a little to understand the contact model and contact stiffness used in YADE._Is it Linear Contact Model at present_? I am trying to explain my question more details as follows,first explanation then my question.

Hello Mohammad.

A linear contact model is anything like F = k*U, provided that k is a constant. All equations you mentionned only differ (or not) in the definition of k, so _they are all linked with a linear contact model._ I'll try to explain the very few ideas on which most equations below are based.

Starting from ref. (03) :

Kn=(kn(A)*kn(B))/(kn(A)+kn(B)))   (eq. 1)

You see here that the stiffness of the contact is in fact the equivalent stiffness of two springs A and B, where ka and kb are the stiffnesses assigned to each individual grains.

Now the question is : how to define ka and kb?
It is in fact convenient to define them as ki = E*ri, with E a constant. The reason is, if you do that, the stiffness of the packing will be independant on the size of particles (you can proove that with dimensional analysis).

Which writes, based on  eq. 1 :

kn = (E.ra*E.rb)/(E.ra+E.rb) = E*ra*rb / (ra+rb)

This is the equation you have in ref. (02) (with "E" named "K", and ra/rb replaced by the mean radius). If you assigned different values of E to different spheres (a mix of hard and soft spheres), it gives the equation in ref. (04) (except the factor 2 that I will comment later, the key point here is just the proportionnality between k and the scaling factor r).

ks, is usually defined as a fraction of kn. Say, for grain i, ks(i)=Vi*kn(i).

From equation 1 again, the equivalent tangential stiffness of the contact will be :

ks = ksa*ksb/(ksa+ksb) , which is, in a condensed form, the equation for ks in ref. (04).

About factor "2" in ref (04) : I added it in the equation so that the Young Modulus Ey of the packing was very close to E in one special case, but this is not something general. Denser/looser packings will have a different proportionality between E and Ey.

If you really want a constant proportionality, you need to replace "2" by a function taking into account the properties of the microstructure. One method for that, based on the coordination number, has been partly implemented in Yade. I can't comment in depth because I did not work on this method, but I know it gave this code with coefficients alpha, beta, gamma, that you can see in MacroMicroElasticRelationship (and not in ref(05) where this is just commented junk code).
So, just take ref (06) as a more sophisticated form of ref (04).

The only expression which - apparently - doesn't fit in this framework, is the one in ref (01). I probably miss some details to understand why. This equation is not used in Yade. But still, the equation in ref (01) defines a constant, so it is relevant for a "linear contact model".

There is nothing fundamental, in fact, in the definition of k. If you want k=10000 N/m at all contacts in your simulations, just write a ExtremellySimpleElasticRelationship engine that will take values of kn and ks as parameters and will assign them to all contacts. There would be nothing wrong in doing that.

Last remark, to prevent some confusion :
I just removed this comment in SimpleElasticRelationships.cpp, because it was just wrong (it was typed by some guys never using PFC apparently ;-)) :

//This is the formula used in PFC-3D
  //Real Kn = 4 * ((Ea+Eb)*0.5) * ((Da+Db)*0.5);
  //Real Ks = Kn/2.0;

I hope it helps.


01) Belheine.N et.al "Numerical simulation of drained triaxial test using 3D discrete element modeling"Computers and Geotechnics,2008 Kn=(kn(A)*kn(B))/(kn(A)+kn(B)))/r page number 03 (Article in press) 02) Kozicki.J et.al."A new open-source software developed for numerical simulations using discrete modeling methods" kn=r*(Kn(A)*Kn(B))/(Kn(A)+Kn(B))) page number 04 03) PFC3D,Theory and Background,Page Number 2-2,Article 2.1.1 Kn=(kn(A)*kn(B))/(kn(A)+kn(B))) kn, normal secant stiffness=Kn,normal tangent stiffness >>My question: what is the significance of "r" average radius,what are the difference among reference 01,02,03. Now inside of YADE (12.rc1) 04) _pkg>dem>engine>engineunit>SimpleElasticRelationships.cpp_ //Real Eab = 2*Ea*Eb/(Ea+Eb);
          //Real Vab      = 2*Va*Vb/(Va+Vb);
Real Dinit = Da+Db; //Real Sinit = Mathr::PI * std::pow( std::min(Da,Db) , 2);

Real Kn = 2*Ea*Da*Eb*Db/(Ea*Da+Eb*Db);//harmonic average of two stiffnesses Real Ks = 2*Ea*Da*Va*Eb*Db*Vb/(Ea*Da*Va+Eb*Db*Va);//harmonic average of two stiffnesses with ks=V*kn for each sphere
// Real alpha = sdecContactModel->alpha;
 //     Real beta       = sdecContactModel->beta;
//      Real gamma      = sdecContactModel->gamma;
// // Real Kn = abs((Eab*Sinit/Dinit)*( (1+alpha)/(beta*(1+Vab) + gamma*(1-alpha*Dab) ) ));
//      Real Ks         = abs(Kn*(1-alpha*Vab)/(1+Vab));
06>pkg>dem>engine>engineunit>MacroMicroElasticRelationships.cpp >>In 04 and 05 the equation is different,can anybody explain a little the difference and significance of 04 and 05. Thanks in advance. Regards,
Mohammad Nurul Islam

Do You Yahoo!?
Tired of spam? Yahoo! Mail has the best spam protection around

Mailing list: https://launchpad.net/~yade-users
Post to     : yade-users@xxxxxxxxxxxxxxxxxxx
Unsubscribe : https://launchpad.net/~yade-users
More help   : https://help.launchpad.net/ListHelp


Chareyre Bruno
Maitre de conference

Grenoble INP
Laboratoire 3SR - bureau E145
BP 53 - 38041, Grenoble cedex 9 - France
Tél : 33 4 56 52 86 21
Fax : 33 4 76 82 70 43