How to impose an initial angular velocity on a moving/free particle?

Asked by kaiset on 2019-06-01

I've tried several ways to impose an initial angular velocity on a free polyhedron but they do not work well as commented in the script and below. Please help me with this. Thanks!
Here are things I've observed: [In the code below, I simulate 2 polyhedrons, poly1 is free and has an initial velocity while poly2 is fixed.]
If I directly use state.angVel or utils.setBodyAngularVelocity to impost an angular velocity on both polyhedrons, nothing happens to poly1, I could observe no rotation but just transnational move. However, I could observe rotation on poly2, which is fixed.
My question is that I want to apply an initial angular velocity on a free polyhedron and see how the initial angular velocity affects its collision with another particle.
If I just fix a polyhedron, and apply both the linear and angular velocity, then I could not observe how it collides with another particle and how the post-impact behavior look like since the fixed one will simply penetrate through another particle, which is not what I want.

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# ENGINES
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Polyhedra_Aabb()],verletDist=0.01),
 InteractionLoop(
  [Ig2_Polyhedra_Polyhedra_PolyhedraGeom()],
  [Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys()],
  [Law2_PolyhedraGeom_PolyhedraPhys_Volumetric()]
 ),
 NewtonIntegrator(damping=0.0,exactAsphericalRot=True,gravity=[0,0,0]),
]

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# POLYHEDRAL PARTICLES
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

from yade import polyhedra_utils
m = PolyhedraMat()
m.density = 2000 #kg/m^3
m.young = 150e6 #Pa
m.poisson = .4
m.frictionAngle = 3.0 #rad

edge=0.025
vertices=[(-edge, -edge, -edge),
          (-edge, edge, -edge),
   ( edge, -edge, -edge),
   ( edge, edge, -edge),
   (-edge, -edge, edge),
          (-edge, edge, edge),
   ( edge, -edge, edge),
   ( edge, edge, edge)]

# Free
b1 = polyhedra_utils.polyhedra(m,v=vertices)
b1.state.pos = (0,0,0)
O.bodies.append(b1)

# Fixed
b2 = polyhedra_utils.polyhedra(m,v=vertices,fixed=True)
b2.state.pos = (4.*edge,0,0)
O.bodies.append(b2)

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# INITIAL VELOCITY
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# Initial translational velocity
O.bodies[0].state.vel=Vector3(1,0,0)

# Initial angular velocity
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
##Directly assigning the angular velocity- Not imposing initial angular velocity- No rotation for the free one but the fixed one will rotate
#O.bodies[0].state.angVel=Vector3(1,1,0)
#O.bodies[1].state.angVel=Vector3(1,1,0)

## using: "utils.setBodyAngularVelocity" - Not working for non-fixed bodies(poly1) but working for poly2
utils.setBodyAngularVelocity(0,5*Vector3(1,1,1))
utils.setBodyAngularVelocity(1,5*Vector3(1,1,1))
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

from yade import qt
v=qt.View()

O.dt=1.0e-8

O.saveTmp()

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
2019-06-10
Last query:
2019-06-10
Last reply:
2019-06-10
Robert Caulk (rcaulk) said : #1

Hello,

Could you be a bit more specific about why these methods "do not work"? Are you seeing an error? If so what is the error? If not what is it you are observing?
I don't see O.run() anywhere in your script, so I don't expect anything to happen in your simulation.
Why are you choosing a manual time step? It is possible the choosen step is too small to observe any significant movements.

Cheers,

Robert

Jan Stránský (honzik) said : #2

Hello,

> ## using: "b.state.angMom" - This could work, but how can I get a specific magnitude of velocity?
> #O.bodies[0].state.angMom=Vector3(b.state.inertia[0]*angVelVector[0],1,1)

have a look at the code how angVel is computed form angMom [1]. After a quick glance, I am not sure if there is an easy inverse of the function..

If there is not, you can state is as an inverse problem ("what andMom is neede to get this angVel") and solve it using some kind of optimization, or at least "brutal force" should work.

cheers
Jan

kaiset (kaiset) said : #3

I have updated the question description.

Best Jan Stránský (honzik) said : #4

Hello,

if you print b.state.angVel, it is set properly. The problem is somewhere in NewtonIntegrator not applying angular motion for fixed aspherical bodies [1], where angMom is is taken into account, but not angVel [1]
Feel free to open a bug report.

A **workaround** is to set
b1.aspherical = False
than the polyhedron does rotate.
I am not sure with the consequences, but **I think** it should be OK for fixed bodies..

cheers
Jan

[1] https://gitlab.com/yade-dev/trunk/blob/master/pkg/dem/NewtonIntegrator.cpp#L242

kaiset (kaiset) said : #5

Thanks Jan Stránský, that solved my question.

kaiset (kaiset) said : #6

Still a **workaround** is to use the angMom

## using: "b.state.angMom" -
Mom=b1.state.ori*b1.state.inertia # This results to angVel=Vector3(1,1,1).
# Modify accordingly, to modify a specific angular velocity. e.g., to get an initial angVel=Vector3(5,6,7)
Mom[0]=Mom[0]*5.
Mom[1]=Mom[1]*6.
Mom[2]=Mom[2]*7.

O.bodies[0].state.angMom=Mom

Janek Kozicki (cosurgi) said : #7

Hi, I moved your bug report to gitlab: https://gitlab.com/yade-dev/trunk/issues/84