Question #698572 on Yade changed:
https://answers.launchpad.net/yade/+question/698572
Rohit John posted a new comment:
Dear Jan,
I performed a simple test on your code. I tried simulating what would
happen a sphere hits the clump "pendulum". Then the angular velocity
about the origin was plotted against time. The system consists of the
clump and the sphere hitting the clump. Since, the external forces are
supposed to act on the hinge, the angular velocity about the hinge
should be conserved. But I see it is not conserved during the collision.
But it is conserved when the clump pendulum is rotating, suggesting your
idea works in that time. It could also be that, there is an error during
the integration of angular momentum. I have pasted the code below.
Kind regards,
Rohit K. John
#
--------------------------------------------------------------------------------------------------------------------------------------------------------------------
Python
from yade import plot
#
----------------------------------------------------------------------------------
bodies
sph1 = sphere((0, -2, 2), .5)
sph2 = sphere((0, 0, 0), .5)
sph3 = sphere((0, 0, 2), .5)
sph3_id = O.bodies.append(sph3)
O.bodies[sph3_id].state.vel = [0,-1,0]
# ---------------------------------------------------------------- clump
clump_id,(sph1_id,sph2_id) = O.bodies.appendClumped((sph1,sph2))
clump = O.bodies[clump_id]
center,radius = Vector3(0,0,0), sqrt(2)
clump.state.refOri = clump.state.ori # not done automatically
clump.state.blockedDOFs = "XYZ" # to force "simple" angular integration
#
----------------------------------------------------------------------------------
engines
O.engines = [
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
),
PyRunner(command = "eccentricAxisEngine()", iterPeriod = 1),
NewtonIntegrator(gravity=(0,0,0), damping=0, label="newton"),
PyRunner(command = "plotter()", iterPeriod = 100)
]
O.dt = .5e-2 * PWaveTimeStep()
#
----------------------------------------------------------------------------------
additional engines
# ----------------------------------------------------------------
constraintClosesPoint
def constraintClosesPoint(pos):
delta = pos - center
direction = delta.normalized()
pos = center + direction * radius
angle = atan2(direction[2],direction[1])
ori = Quaternion()
ori.setFromTwoVectors(Vector3.UnitX,direction)
return pos,ori
# ----------------------------------------------------------------
eccentricAxisEngine
def eccentricAxisEngine():
dt = O.dt
s = clump.state
f = O.forces.f(clump_id)
a = f / s.mass + newton.gravity
pos1 = s.pos + (s.vel+a*dt)*dt
pos0,ori0 = constraintClosesPoint(pos1)
s.vel = (pos0 - s.pos) / dt - a*dt
#
v1 = s.ori * Vector3.UnitX
v2 = ori0 * Vector3.UnitX
ang = v1.cross(v2)
ang = ang / O.dt
s.angVel = ang
# ---------------------------------------------------------------- plotter
plot.plots = {'t':('ax', 'ay', 'az')}
def plotter():
body_ids = [clump_id, sph3_id]
totalangVel = Vector3([0,0,0])
for i in body_ids:
totalangVel = totalangVel + getAngVelAboutOrigin(i)
plot.addData(t = O.time, ax = totalangVel[0],
ay = totalangVel[1], az = totalangVel[2]
)
plot.plot()
# ----------------------------------------------------------------
getAngVelAboutOrigin
def getAngVelAboutOrigin(id):
state = O.bodies[id].state
pos = state.pos
vel = state.vel
angMom = state.angMom
mass = state.mass
total_ang_mom = angMom + mass*pos.cross(vel)
return total_ang_mom
--
You received this question notification because your team yade-users is
an answer contact for Yade.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
