yade-users team mailing list archive

[Rohit John <question698572@xxxxxxxxxxxxxxxxxxxxx>]

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.