i.geom.normal is not normalized (unit) vector from the branch vector (for the case of sphere-sphere contact)

Asked by Zhicheng Gao

I want to calculate the mechanical anisotropy, the branch vector is needed. However, I find the contact normal obtained by i.geom.normal is not normalized (unit) vector from the branch vector (for the case of sphere-sphere contact), The specific process is as follows:
(1) load saved file by "O.load('compactedStatesuffusion25.yade.gz')" ;
(2) find all the interactions of particle 12 by "for i in O.bodies[12].intrs(): print (i.isReal,i.id1,i.id2)" and the result is as follows:
(True, 3, 12)
(True, 12, 97)
(True, 12, 211)
(True, 427, 12)
(True, 12, 1033)
(True, 1171, 12)
(True, 1196, 12)
(True, 12, 2736)
(True, 2777, 12)
(True, 12, 3222)
(True, 3453, 12)
(True, 4202, 12)
(True, 12, 7845)
(True, 12, 8119)
(True, 8496, 12)
(True, 10817, 12)
(True, 11570, 12)
(True, 12, 13485)
(True, 12, 15013)
(True, 12, 19529)
(True, 12, 21199)
(True, 12, 22291)
(True, 12, 22449)
(True, 12, 22803)
(True, 12, 23442)
(True, 12, 24412)
(True, 25196, 12)
(True, 12, 29244)
(True, 12, 30387)
(True, 12, 34784)
(True, 12, 35280)
(True, 36427, 12)
(True, 37206, 12)
(True, 12, 40297)
(True, 41898, 12)
(True, 12, 42413)
(True, 12, 44052)
(True, 12, 45495)
(True, 45948, 12)
(True, 12, 47578)
(True, 12, 48322)
(True, 12, 49112)
(True, 12, 49738)
(True, 49923, 12);
(3) then print the branch vector, contact normal vector and contact normal force vector of O.interactions[3,12] by the follow code:
i = O.interactions[3,12]
b1,b2 = [O.bodies[j] for j in (i.id1,i.id2)]
print("branch12",b2.state.pos - b1.state.pos)
print("branch21",b1.state.pos - b2.state.pos)
print("normal",i.geom.normal)
print("force",i.phys.normalForce + i.phys.shearForce)
print("force1",O.forces.f(i.id1))
print("force2",O.forces.f(i.id2))

the result is as follows:

('branch12', Vector3(0.005264805042519041,-0.0011961713201094834,-0.0024246928681093648))
('branch21', Vector3(-0.005264805042519041,0.0011961713201094834,0.0024246928681093648))
('normal', Vector3(-0,-1,-0))
('force', Vector3(0,-0.6724835758663493,-0))
('force1', Vector3(0,28.188489135826057,0))
('force2', Vector3(-0.0003050240398732389,0.00005915155744011741,0.00009299475981276988))

It can be seen that the branch vector, contact normal vector and contact normal force vector are not in the same direction. Why is this?
Thanks!

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Zhicheng Gao
Solved:
Last query:
Last reply:
Revision history for this message
Jan Stránský (honzik) said :
#1

Hello,

there might be several reasons, difficult to tell without actual data [1]. For a general case:

A) contact normal vector vs. contact normal force:
> ... contact normal force vector ...
> print("force",i.phys.normalForce + i.phys.shearForce)

you are printing TOTAL force, not normal force.
If i.geom.shearForce is non-zero, then you naturally get non-parallel vector

O.forces.f(i.id1) result is "independent" of constituent interactions

B) branch vector vs. contact normal vector:

i.geom.normal is computed in InteractionLoop.
Body positions are THEN updated by NewtonIntegrator.
So b2.state.pos - b1.state.pos and i.geom.normal (in general) are not parallel.

Cheers
Jan

[1] https://www.yade-dem.org/wiki/Howtoask

Revision history for this message
Zhicheng Gao (zhichenggao) said :
#2

Thanks, Jan

>you are printing TOTAL force, not normal force. If i.geom.shearForce is non-zero, then you naturally get non-parallel vector
I have modified it to print("force",i.phys.normalForce)and I found out that it was my mistake, body 3 is a wall, not a particle. I checked the other contacts and they all met the parallel conditions. The question now is why this phenomenon occurs on the wall.

Best wishes
Gao

Revision history for this message
Jan Stránský (honzik) said :
#3

> The question now is why this phenomenon occurs on the wall.

just print the wall's state.pos and it should answer you itself :-)
Basically, wall's state.pos has nothing to do with the contact point, which determines geom.normal

cheers
Jan

Revision history for this message
Zhicheng Gao (zhichenggao) said :
#4

For i = O.interactions[A,B],is i.phys.normalForce the force exerted on grain A by grain B or grain B by grain A?

Revision history for this message
Jan Stránský (honzik) said :
#5

a) see documentation [2]
b) have a try yourself
Cheers
Jan

[2] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.FrictPhys.normalForce

Revision history for this message
Zhicheng Gao (zhichenggao) said :
#6

Jan, thank you very much for your help and quick reply!