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

Best

Vincent

## Question information

Language:
English Edit question
Status:
Solved
For:
Assignee:
No assignee Edit question
Solved by:
Rioual
Solved:
2020-04-24
Last query:
2020-04-24
2020-04-21
 Jan Stránský (honzik) said on 2020-04-20: #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

 Jan Stránský (honzik) said on 2020-04-20: #2
 Rioual (francois-rioual-v) said on 2020-04-21: #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:

How can I get it for my purpose ??

Best

Vincent

 Jan Stránský (honzik) said on 2020-04-21: #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:

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 (francois-rioual-v) said on 2020-04-21: #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 on 2020-04-21: #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.shearIncrement... are just one vector.
If you need/want to decompose it into two in-tangential-plane 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 (francois-rioual-v) said on 2020-04-21: #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 on 2020-04-21: #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

cheers
Jan

 Rioual (francois-rioual-v) said on 2020-04-21: #9

Hello Jan,

To be more precise:

F_n = f(alfa) * normal_strain_rate + g(alfa)

F_t = h(alfa) * tangential_strain_rate

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 on 2020-04-21: #10

> F_n = f(alfa) * normal_strain_rate + g(alfa)
> F_t = h(alfa) * tangential_strain_rate

# 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_strain_rate = i.geom.shearIncrement / O.dt / length
F_n = f(alfa) * normal_strain_rate + g(alfa)
F_t = h(alfa) * tangential_strain_rate
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 (francois-rioual-v) said on 2020-04-22: #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
Maybe treated elsewhere?)

V.

 Jan Stránský (honzik) said on 2020-04-22: #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:
- 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 (francois-rioual-v) said on 2020-04-22: #13

Thanks Jan, that was great !
best wishes
V.

 Rioual (francois-rioual-v) said on 2020-04-24: #14

Thanks for all,

V.