yade-users team mailing list archive

[Rohit John <question698572@xxxxxxxxxxxxxxxxxxxxx>]

Question #698572 on Yade changed:
https://answers.launchpad.net/yade/+question/698572

    Status: Answered => Open

Rohit John is still having a problem:
Dear Jan,

I tried the same test, explain in the previous comment, on a clump which
is force to rotate about the origin using a penalty force model. The
change in angular momentum is still present, however, the magnitude is
much less. I suppose this could be because of an error in the
integration of the angular momentum.

Kind regards,
Rohit K. John

# -------------------------------------------------------------------------------------------------------------------------------------------------------------------- Python
from math import *
from yade import *
from yade import utils, geom, plot


# ---------------------------------------------------------------------------------- bodies
sph1 = sphere((0, -2, 2), .5)
sph2 = sphere(center=(0, 0, 0), radius=.5, fixed=True)
sph3 = sphere((0, 0, 2), .5)

sph1_id = O.bodies.append(sph1)
sph2_id = O.bodies.append(sph2)

sph3_id = O.bodies.append(sph3)
O.bodies[sph3_id].state.vel = [0,-1,0]

g        =  Vector3([0,0,0])
stiffnes = 1e7
# ---------------------------------------------------------------- clump
clump_id = O.bodies.clump([sph1_id, sph2_id])


# ---------------------------------------------------------------------------------- engines
O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb()]),
        InteractionLoop(
                [Ig2_Sphere_Sphere_ScGeom()],  # collision geometry
                [Ip2_FrictMat_FrictMat_FrictPhys()],  # collision "physics"
                [Law2_ScGeom_FrictPhys_CundallStrack()]  # contact law -- apply forces
        ),
        PyRunner(command = "eccentricAxisEngine()", iterPeriod = 1),
        NewtonIntegrator(gravity=g, damping=0.0),
        PyRunner(command = "plotter()", iterPeriod = 100)
]

#
----------------------------------------------------------------------------------
additional engines

def eccentricAxisEngine():
    pos_0 = O.bodies[sph2_id].state.pos
    force = -stiffnes*pos_0
    O.forces.addF(id = sph2_id, f = force)

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()

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
# ---------------------------------------------------------------------------------- sim setup
O.dt = .5e-2 * PWaveTimeStep()
O.saveTmp()

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.