yade-users team mailing list archive
-
yade-users team
-
Mailing list archive
-
Message #01266
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.
Bruno
Reference:
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
05)_pkg>dem>engine>standaloneengine>GlobalStiffnessTimestepper.cpp_
// 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
http://mail.yahoo.com
------------------------------------------------------------------------
_______________________________________________
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
________________
References
