Using O.forces along normal and tangential directions at contact
Hello,
I would like to add a force at the contact particles with O.forces (doing a loop on all interactions between particles).
This force evolves in a complex way in time and has a normal component Fn dependant on normal strain rate and tangential component Ft dependant on tangential strain rate. Can I have access easily normal strain rates (U_n) , tangential strain rates ( U_t) and components of the normal and tangential unit vectors ( n and t) in global coordinates (to use O.forces) ??? by ScGeom ? How ??
Thanks for your help,
Best
Vincent
Question information
 Language:
 English Edit question
 Status:
 Solved
 For:
 Yade Edit question
 Assignee:
 No assignee Edit question
 Solved by:
 Rioual
 Solved:
 20200424
 Last query:
 20200424
 Last reply:
 20200421
Jan Stránský (honzik) said :  #1 
Hello,
> normal strain rates (U_n)
geom.incidentVel [1], together with geom.refR1[2] and geom.refR2 [3] to convert velocity to strain rate
> tangential strain rates ( U_t)
geom.shearInc [4], together with O.dt to convert strain increment to strain rate
> components of the normal and tangential unit vectors ( n and t) in global coordinates
n = geom.normal [5]
t is more tricky, there is no unique "t unit vectors", but only the tangential plane is defined uniquely. You can choose some two perpendicular vectors arbitrarily if you really need it (often working with global coordinate system but assuming vectors being only in the tangential plane is suitable).
cheers
Jan
[1] https:/
[2] https:/
[3] https:/
[4] https:/
[5] https:/
Jan Stránský (honzik) said :  #2 
Rioual (francoisrioualv) said :  #3 
Dear Jan,
Thanks for your precious input but i do not get you concerning the determination of t
being tricky.
As far as I know, I thought it was precisely defined in DEM in order to compute contact forces:
https:/
How can I get it for my purpose ??
Best
Vincent
Jan Stránský (honzik) said :  #4 
> but i do not get you concerning the determination of t being tricky.
Maybe "tricky" is not the right word..
the point is that for the contact, you have geom.normal. So the "n" is uniquely defined. Then you have a plane perpendicular to n, which is uniquely defined. Then in that plane, you can choose arbitrarily two perpendicular vectors to be "t1" and "t2" (out of infinitely many possibilities).
This is the "tricky" part.
IMO if you do not really really need it (which should be the case) or do not know exactly what you are doing (you would not ask), then do not use shear components in local coordinate system.
> As far as I know, I thought it was precisely defined in DEM in order to compute contact forces:
> https:/
contact forces are simply computed in global coordinate system. This is exactly an example of my remark:
>> (often working with global coordinate system but assuming vectors being only in the tangential plane is suitable).
the picture is a 2D illustrative representation (in 2D the situation is much easier), the computations below it are for the general 3D case.
> How can I get it for my purpose ??
if your purpose is adding forces with O.forces, then use the "global coordinate system approach" (as O.forces are in global coordinates anyway).
Or
>> You can choose some two perpendicular vectors arbitrarily if you really need it
cheers
Jan
Rioual (francoisrioualv) said :  #5 
Dear Jan,
Thanks, i get it but how do i practically
do to get the components of t in global
coordinates ??
That's the remaining bit for my problem..
Sorry for not being clear,
Best
Vincent
Jan Stránský (honzik) said :  #6 
> IMO if you do not really really need it (which should be the case) or do not know exactly what you are doing (you would not ask), then do not use shear components in local coordinate system.
Also sorry for not being clear, the best is (IMO) not to use shear components at all and simply work in global coordinate system. E.g. the phys.shearForce, geom.shearIncre
If you need/want to decompose it into two intangentialplane components, it is possible, but, as mentioned above, often (always?) it is not necessary.
If it is still not clear, maybe you can describe more in detail "evolves in a complex way in time" for "if you do not really really need it" verification.
cheers
Jan
Rioual (francoisrioualv) said :  #7 
Hello Jan,
I have to implement a supplementary force at each contact which can defined as:
F_n = function1 of (normal strain rate, ...)
F_t = function2 of (tangential strain rate, ..)
where normal direction is along the centers of the particles and tangential direction is defined according
to the relative velocity of the particles at the contact....
Now, If I want to use O.forces, I need to convert these components F_n and F_t in global coordinates.
I don't see how I can avoid that task thanks to a vectorial relationship in my case !!
...Just wanted to know if there was an easy way (functions) to compute that in yade..
??
Best wishes !
V.
Jan Stránský (honzik) said :  #8 
> F_n = function1 of (normal strain rate, ...)
> F_t = function2 of (tangential strain rate, ..)
this is already clear from from the original post :) I think the discussion needs more information on this.
Is it possible to provide more details, like a specific formula, on what other quantities the result depend etc.?
> and tangential direction is defined according to the relative velocity of the particles at the contact....
What is "the relative velocity of the particles at the contact"? something else then geom.incidentVel and geom.shearIncrement (or their strain rate or velocity counterparts)?
> Now, If I want to use O.forces, I need to convert these components F_n and F_t in global coordinates.
Depending on the meaning of "components F_n and F_t" (??):
1) normal + shear components, simply:
F = F_n + F_t
2) shear components
you don't (shouldn't) need to convert anything. You have all variables as vectors in global coordinate system, O.forces are vectors in global coordinate system, so just compute with vectors in global coordinate system (that's why I require more info on actual computation, I don't see here any use case where shear components would be needed)
> I don't see how I can avoid that task thanks to a vectorial relationship in my case !!
The vectorial relationshop is actually the way to avoid it..
> ...Just wanted to know if there was an easy way (functions) to compute that in yade.. ??
see above
waiting for more info
cheers
Jan
Rioual (francoisrioualv) said :  #9 
Hello Jan,
To be more precise:
F_n = f(alfa) * normal_strain_rate + g(alfa)
F_t = h(alfa) * tangential_
alfa is a dynamic variable which varies at each time step according to an evolution differential equation
that involves the geometry of contact and normal strain rate...etc
(d(alfa)/dt = function(alfa, normal_strain_rate, geometrical parameters , physical parameters)....
....So that I want to use O.forces to implement that in an incremental way !!
When i say "relative velocities of the particles at contact", I just said "at contact" because of potential spin on the particles...
All the best!
Vincent
Jan Stránský (honzik) said :  #10 
> F_n = f(alfa) * normal_strain_rate + g(alfa)
> F_t = h(alfa) * tangential_
thanks for more info, then it is as easy as:
# i is an Interaction instance, e.g. within "for i in O.interactions:"
length = i.geom.refR1 + i.geom.refR2 # or something similar
normal_strain_rate = i.geom.incidentVel / length
tangential_
F_n = f(alfa) * normal_strain_rate + g(alfa)
F_t = h(alfa) * tangential_
F = F_n + F_t
# add F (and probably corresponding torque, too) to O.forces for both bodies (maybe open a new question for this topic if needed)
just a note, all this seems like a task for a Law2 functor (yes, which would need to be written in C++ and compiled with Yade, but the computation would be significantly faster....).
So according to your problem, plans, domain size...., consider it as a future step.
cheers
Jan
Rioual (francoisrioualv) said :  #11 
Hello jan,
This becomes very clean now, thanks.
What is the problem in doing addF for
both bodies ? You propose me to open a
New question...There is a short and efficient
way tu do it, so ?
 concerning the suggestion of task for a
Law2 functor, this is definitely a future
step for me as i will have lots of particles
to deal with, can you advice me an easy way
to handle that as i have no expérience of
it and i start with yade ?? (A basic question
Maybe treated elsewhere?)
Thank you for your enlightements,
V.
Jan Stránský (honzik) said :  #12 
Concerning O.forces, there are a few tricky points:
 the force should be opposite on both bodies
 addF/setPermF (depending on Yade version) overwrites existing values, one need to take care if dealing with multiple bodies
 maybe something else..
 same for torque
just to make the discussion clean, please open a new question in case of a new question :)
Concerning Law2  it is not too easy task but not too difficult at the same time). There are a few steps to take:
 compile yade ffom source
 decide if need new Material and/or new IPhys or just new Law2

 it is good to create a small "one interaction" benchmark (or a few of them) with known response and checking the implemented model with the known values
 ...
again, a new question(s) is appropriate when the time comes :)
cheers
Jan
Rioual (francoisrioualv) said :  #13 
Thanks Jan, that was great !
best wishes
V.
Rioual (francoisrioualv) said :  #14 
Thanks for all,
V.