yade-users team mailing list archive
-
yade-users team
-
Mailing list archive
-
Message #24123
[Question #693479]: How to initialize the angular velocity of a box made of pfacets
New question #693479 on Yade:
https://answers.launchpad.net/yade/+question/693479
Hello,
I created a box using pFacets. Now I want to initialise its angular velocity. I did so by assigning the same angular velocity to all bodies in the pfacet box(nodes, cylinders and pfacet) using O.bodies[i].state.angVel, and I also assigned the appropriate linear velocities using v = r x w, where v is the velocity, r is the radius vector and w is the angular velocity. Gravity is 0 and there is no other interaction in the simulation. I ran the code and saw the box rotating in the viewport without any erratic motion. However, it is behaving differently when I plot angVel for a body in the pfacet box. It was not constant. I even tried averaging this value for all bodies in the pfacet box. It did not work. The code for this case is given below as Code A.
To check if the angVel property has an issue I tried the following. I set two spheres, positioned close to two opposite edges and moving towards the box in opposite direction, to collide with the box. This did the trick and made the box rotate. When I plotted angVel it shows a constant (almost) value. The code for this case is given below as Code B.
My question is there a way to set the angular velocity. I do not want to use angMom as the moment of inertia for all bodies in the pfacet box are different and I guess some objects are planar so they have 0 moments of inertia.
############################################################################################################## Code A
from yade import *
from yade import plot
from yade.gridpfacet import *
pfacet_rad = 0.005
young = 1e6
density = 1e3
frictAng = radians(30)
poisson = 0.3
# function --------------------------------------------------------------------------------------
## box creator function -----------------------------------------------------------
def symmetric_cuboid_pfacet(CM, length, breadth, height, radius, int_mat, ext_mat, color = [0.5, 0.5, 0.5]):
'''
makes a symmetric cuboid using pfacet.
'''
l = length
b = breadth
h = height
cm1, cm2, cm3 = CM # the centre of mass of the cube
# coordinates of the vertices
vertices = [
[cm1 + l/2, cm2 + b/2, cm3 + h/2],
[cm1 + l/2, cm2 - b/2, cm3 + h/2],
[cm1 - l/2, cm2 - b/2, cm3 + h/2],
[cm1 - l/2, cm2 + b/2, cm3 + h/2],
[cm1 + l/2, cm2 + b/2, cm3 - h/2],
[cm1 + l/2, cm2 - b/2, cm3 - h/2],
[cm1 - l/2, cm2 - b/2, cm3 - h/2],
[cm1 - l/2, cm2 + b/2, cm3 - h/2],
[cm1 + l/2, cm2 , cm3 ],
[cm1 - l/2, cm2 , cm3 ],
[cm1 , cm2 + b/2, cm3 ],
[cm1 , cm2 - b/2, cm3 ],
[cm1 , cm2 , cm3 + h/2],
[cm1 , cm2 , cm3 - h/2],
]
# connecting the vertices to form faces
face = [
[13, 2, 1], [13, 3, 2], [13, 4, 3], [13, 1, 4],
[14, 5, 6], [14, 6, 7], [14, 7, 8], [14, 8, 5],
[11, 1, 5], [11, 5, 8], [11, 8, 4], [11, 4, 1],
[12, 2, 3], [12, 3, 7], [12, 7, 6], [12, 6, 2],
[9, 1, 2], [9, 2, 6], [9, 6, 5], [9, 5, 1],
[10, 3, 4], [10, 4, 8], [10, 8, 7], [10, 7, 3],
]
# generating the grid stuff
nodesIds = []
cylIds = []
pFacet = []
for i in face:
pfacetCreator1(
[
vertices[i[0] - 1],
vertices[i[1] - 1],
vertices[i[2] - 1]
],
radius,
nodesIds = nodesIds,
cylIds = cylIds,
pfIds = pFacet,
wire = False,
fixed = False,
color = color,
materialNodes = int_mat,
material = ext_mat,
)
return nodesIds ,cylIds, pFacet
## cross product function ---------------------------------------------------------
def cross_prd(a,b):
'''
Returns the cross product of the vectors a and b
'''
a1,a2,a3 = a
b1,b2,b3 = b
c1 = a2*b3 - a3*b2
c2 = a3*b1 - a1*b3
c3 = a1*b2 - a2*b1
c = [c1 ,c2, c3]
return c
## vector difference function -----------------------------------------------------
def subtract(a,b):
'''
returns the difference of the vectors a and b
'''
a1,a2,a3 = a
b1,b2,b3 = b
c = [a1-b1, a2-b2, a3-b3]
return c
## add an angular velocity to a list of bodies function ---------------------------
def body_angular_velocity_add(ang_vel, body_list, CM):
'''
Adds a value to the the initial angular velocity (ang_vel) of the pfacet body (body_list) about
the centre point CM
'''
for i in body_list:
pos = O.bodies[i].state.pos
r = subtract(pos, CM)
vel = cross_prd(r, ang_vel)
O.bodies[i].state.vel = O.bodies[i].state.vel + vel
O.bodies[i].state.angVel = O.bodies[i].state.angVel + ang_vel
# material ------------------------------------------------------------------------------------------
## pfacet internal material ------------------------------------------------------
O.materials.append(
FrictMat(
young = young,
poisson = poisson,
density = density,
frictionAngle = frictAng,
label = 'pFacet_out',
)
)
## pfacet external material ------------------------------------------------------
O.materials.append(
CohFrictMat(
young = young,
poisson = poisson,
density = density,
frictionAngle = frictAng,
momentRotationLaw = True,
normalCohesion = 1e40,
shearCohesion = 1e40,
label = 'pFacet_in'
)
)
# engines -------------------------------------------------------------------------------------------
O.engines = [
ForceResetter(),
InsertionSortCollider([
Bo1_GridConnection_Aabb(),
Bo1_PFacet_Aabb()
]),
InteractionLoop(
[
Ig2_GridConnection_GridConnection_GridCoGridCoGeom(),
Ig2_GridNode_GridNode_GridNodeGeom6D(),
Ig2_GridConnection_PFacet_ScGeom(),
Ig2_PFacet_PFacet_ScGeom()
],
[
Ip2_FrictMat_FrictMat_FrictPhys(),
Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow = True, setCohesionOnNewContacts = False)
],
[
Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(),
Law2_ScGeom_FrictPhys_CundallStrack(),
Law2_ScGeom6D_CohFrictPhys_CohesionMoment()
],
),
NewtonIntegrator(gravity = (0,0,0),damping = 0.0),
PyRunner(command='addPlot()',iterPeriod=10)
]
O.dt = 1e-6
# bodies --------------------------------------------------------------------------------------------
## box ----------------------------------------------------------------------------
(
node_ids,
cyl_ids,
pFacet_ids
) = symmetric_cuboid_pfacet(
[0,0,0],
length =0.05,
breadth =0.05,
height =0.05,
radius =pfacet_rad,
int_mat = 'pFacet_in',
ext_mat = 'pFacet_out',
color = [0.5, 0.5, 0.5]
)
## defining the velocity ----------------------------------------------------------
body_angular_velocity_add(
[10,0,0],
node_ids + cyl_ids + pFacet_ids,
[0,0,0]
)
# plots ---------------------------------------------------------------------------------------------
plot.plots = {'t':'wx'}
def angvelx():
'''Returns the average angular velocity along the x axis'''
wx = 0
for i in (node_ids + cyl_ids + pFacet_ids):
wx = wx + O.bodies[i].state.angVel[0]
n = len(node_ids + cyl_ids + pFacet_ids)
return wx/n
def addPlot():
plot.addData(
t = O.time,
wx = angvelx(),
)
plot.plot()
############################################################################################################## End of code A
############################################################################################################## Code B
# The code for box rotated by impacting it with spheres
# Code B
# new added plots
from yade import *
from yade import plot
from yade.gridpfacet import *
# set up sims ---------------------------------------------------------------------------------------
pfacet_rad = 0.005
young = 1e9
density = 1e3
frictAng = radians(30)
poisson = 0.3
# function --------------------------------------------------------------------------------------
## box creator function -----------------------------------------------------------
def symmetric_cuboid_pfacet(CM, length, breadth, height, radius, int_mat, ext_mat, color = [0.5, 0.5, 0.5]):
'''
makes a symmetric cuboid using pfacet.
'''
l = length
b = breadth
h = height
cm1, cm2, cm3 = CM # the centre of mass of the cube
# coordinates of the vertices
vertices = [
[cm1 + l/2, cm2 + b/2, cm3 + h/2],
[cm1 + l/2, cm2 - b/2, cm3 + h/2],
[cm1 - l/2, cm2 - b/2, cm3 + h/2],
[cm1 - l/2, cm2 + b/2, cm3 + h/2],
[cm1 + l/2, cm2 + b/2, cm3 - h/2],
[cm1 + l/2, cm2 - b/2, cm3 - h/2],
[cm1 - l/2, cm2 - b/2, cm3 - h/2],
[cm1 - l/2, cm2 + b/2, cm3 - h/2],
[cm1 + l/2, cm2 , cm3 ],
[cm1 - l/2, cm2 , cm3 ],
[cm1 , cm2 + b/2, cm3 ],
[cm1 , cm2 - b/2, cm3 ],
[cm1 , cm2 , cm3 + h/2],
[cm1 , cm2 , cm3 - h/2],
]
# connecting the vertices to form faces
face = [
[13, 2, 1], [13, 3, 2], [13, 4, 3], [13, 1, 4],
[14, 5, 6], [14, 6, 7], [14, 7, 8], [14, 8, 5],
[11, 1, 5], [11, 5, 8], [11, 8, 4], [11, 4, 1],
[12, 2, 3], [12, 3, 7], [12, 7, 6], [12, 6, 2],
[9, 1, 2], [9, 2, 6], [9, 6, 5], [9, 5, 1],
[10, 3, 4], [10, 4, 8], [10, 8, 7], [10, 7, 3],
]
# generating the grid stuff
nodesIds = []
cylIds = []
pFacet = []
for i in face:
pfacetCreator1(
[
vertices[i[0] - 1],
vertices[i[1] - 1],
vertices[i[2] - 1]
],
radius,
nodesIds = nodesIds,
cylIds = cylIds,
pfIds = pFacet,
wire = False,
fixed = False,
color = color,
materialNodes = int_mat,
material = ext_mat,
)
return nodesIds ,cylIds, pFacet
# material ------------------------------------------------------------------------------------------
## pfacet internal material ------------------------------------------------------
O.materials.append(
FrictMat(
young = young,
poisson = poisson,
density = density,
frictionAngle = frictAng,
label = 'pFacet_out',
)
)
## pfacet external material ------------------------------------------------------
O.materials.append(
CohFrictMat(
young = 1000*young,
poisson = poisson,
density = density,
frictionAngle = frictAng,
momentRotationLaw = True,
normalCohesion = 1e40,
shearCohesion = 1e40,
label = 'pFacet_in'
)
)
## sphere material ---------------------------------------------------------------
O.materials.append(
FrictMat(
young = young,
poisson = poisson,
density = 10000*density,
frictionAngle = frictAng,
label = 'sphere',
)
)
# engines -------------------------------------------------------------------------------------------
O.engines = [
ForceResetter(),
InsertionSortCollider([
Bo1_Sphere_Aabb(),
Bo1_GridConnection_Aabb(),
Bo1_PFacet_Aabb(),
]),
InteractionLoop(
[
Ig2_Sphere_Sphere_ScGeom(),
Ig2_GridConnection_GridConnection_GridCoGridCoGeom(),
Ig2_GridNode_GridNode_GridNodeGeom6D(),
Ig2_GridConnection_PFacet_ScGeom(),
Ig2_PFacet_PFacet_ScGeom(),
Ig2_Sphere_PFacet_ScGridCoGeom(),
Ig2_Sphere_GridConnection_ScGridCoGeom(),
],
[
Ip2_FrictMat_FrictMat_FrictPhys(),
Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow = True, setCohesionOnNewContacts = False)
],
[
Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(),
Law2_ScGeom_FrictPhys_CundallStrack(),
Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),
Law2_ScGridCoGeom_FrictPhys_CundallStrack(),
],
),
NewtonIntegrator(gravity = (0,0,0),damping = 0.0),
PyRunner(command='addPlot()',iterPeriod=10)
]
# bodies --------------------------------------------------------------------------------------------
## box ----------------------------------------------------------------------------
(
node_ids,
cyl_ids,
pFacet_ids
) = symmetric_cuboid_pfacet(
[0,0,0],
length =0.05,
breadth =0.05,
height =0.05,
radius =pfacet_rad,
int_mat = 'pFacet_in',
ext_mat = 'pFacet_out',
color = [0.5, 0.5, 0.5]
)
## box ----------------------------------------------------------------------------
## spheres ------------------------------------------------------------------------
s3 = utils.sphere([0, 0.04, 0.025],0.005,material='sphere')
s4 = utils.sphere([0,-0.04,-0.025],0.005,material='sphere')
s3id = O.bodies.append(s3)
s4id = O.bodies.append(s4)
O.bodies[s3id].state.vel = [0,-1,0]
O.bodies[s4id].state.vel = [0,1,0]
O.saveTmp()
# plots ---------------------------------------------------------------------------------------------
plot.plots = {'t':'wx'}
def angvelx():
'''Returns the average angular velocity along the x axis'''
wx = 0
for i in (node_ids + cyl_ids + pFacet_ids):
wx = wx + O.bodies[i].state.angVel[0]
n = len(node_ids + cyl_ids + pFacet_ids)
return wx/n
def addPlot():
plot.addData(
t = O.time,
wx = angvelx(),
)
plot.plot()
O.dt = utils.PWaveTimeStep()
O.saveTmp()
--
You received this question notification because your team yade-users is
an answer contact for Yade.
