# 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.

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

_______________________________________________
```
```

--

_______________
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
________________

```