Yade

ViscoEl contact normal stiffness kn in the source file is not consistent with that in the documentation

Asked by Yuxuan Wen

Hello,

I am trying to use the visco-elastic model to do the simulation. I checked the kn calculation in the source file and that in the documentation, then I found that there is a factor 2 difference.

In the documentation, the stiffness is kn=2*kn1*kn2/(kn1+kn2), as shown in the following link:
https://yade-dem.org/doc/yade.wrapper.html?highlight=law2_scgeom_viscelphys_basic#yade.wrapper.Law2_ScGeom_ViscElPhys_Basic

However, in the source file, the Kn is calculated by kn= kn1*kn2/(kn1+kn2), as shown in source file line 279 and 315-323:
https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/dem/ViscoelasticPM.cpp
279: phys->kn = contactParameterCalculation(kn1, kn2)

315: Real contactParameterCalculation(const Real& l1, const Real& l2)
316: {
317: // If one of paramaters > 0. we DO NOT return 0
318: Real a = (l1 ? 1 / l1 : 0) + (l2 ? 1 / l2 : 0);
319: if (a)
320: return 1 / a;
321: else
322: return 0;
323: }

Why they are different?

I know in the DEM configuration, Fn=kn1*dl1=kn2*dl2=kn*(dl1+dl2), and then we get kn=kn1*kn2/(kn1+kn2) by derivation , but this has a requirement that (dl1+dl2) is the length corresponds to kn. However, the penetration depth in YADE is calculated by dl=R1+R2-L where L is the length from sphere1 center to sphere2 center, in this case, dl=1/2*(dl1+dl2). Then F=kn*(dl1+dl2)=kn*(2*dl)=2*kn*dl=2*kn1*kn2/(kn1+kn2)*dl, and the new kn = 2*kn1*kn2/(kn1+kn2), that is exactly the kn expression in documentation. However, in the source file, this factor 2 is missing. Am I thinking correctly?

Would you please help me? Thank you!

Kind Rrgards,
Yuxuan Wen

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jérôme Duriez
Solved:
Last query:
Last reply:
Revision history for this message
Launchpad Janitor (janitor) said : 		#1

This question was expired because it remained in the 'Open' state without activity for the last 15 days.

Revision history for this message
Jérôme Duriez (jduriez) said : 		#2

Hi,

If your question was "there is a mismatch between documentation and source code, which should I trust?" the answer is "source code".
Behavior of YADE executable is obviously the one of its source code (the one of the used version). Documentation just attempts to retranscript this behavior in a non-C++, more friendly, language. Sometimes the attempt fails.. (if that's the cas, it may be fixed in a near future)

If your question is "how to justify the kn constant expression, with or without a 2", this may have been discussed few times on the mailing list in the last 10 years (..) and I think there is simply is no answer.
As far as I'm concerned, I can imagine no bullet-proof theoretical development for such a linear approximation (why Fn = kn dl1 would hold ? This is not what Hertz theory says). And there are then several (size-independent, at the macro scale) expressions that are all equally as good, with or without multiplying factors 2, or pi, or whatever number you could think of.

Revision history for this message
Yuxuan Wen (wenyuxuan) said : 		#3

Thank you Jerome!

So do you mean that a constant factor such as 2, or pi, doesn't really affect the results of simulation? Because the fundamental equation in DEM Fn=kn dl1 is also not strictly proved and is just a approximation.

Kind Regards,
Yuxuan

Revision history for this message
Best Jérôme Duriez (jduriez) said : 		#4

Even though having a different kn will affect the simulation, I mean that you can always get the same kn (and the same simulation results), as soon as you change your kn-related user-parameter (that will combine with 2 or pi to give kn).

And there is no physical grounds to discriminate between one user-parameter value (the one combined with 2) or another (the one combined with pi) => no reason to prefer the expression with 2 or the expression with pi.

Revision history for this message
Yuxuan Wen (wenyuxuan) said : 		#5

Thanks Jérôme Duriez, that solved my question.