New question #699510 on Yade:
https://answers.launchpad.net/yade/+question/699510
Hi,
I'd like to ask that how to get the resultant force on a particle.
I came out with two ways to get the (magnitude of) resultant force on a
particle as follow.
(1) using O.forces.f(id)
(2) calculating from the interactions on that particle:
a=[]
for i in O.bodies[id].intrs():
f = i.phys.normalForce + i.phys.shearForce
a.append(f)
resultF=Vector3(0,0,0)
for i in range(len(a)):
resultF=resultF+a[i]
However, in some of situations, these two ways can have the same magnitude of
resultant force, for instance, using MWE1, the magnitude of resultant force
obtained from the two approaches are the same. But using MWE2, I got different
results from the two approaches. (both the two MWE are slightly modified from
YADE example[1])
Could you please help me with: whether the way I calculate the resultant force
is correct; if they are correct, why we can get different results?
########### MWE1 which can lead to the same results###########
O.bodies.append(
[
# fixed: particle's position in space will not change (support)
sphere(center=(0, 0, 0), radius=.5, fixed=True),
# this particles is free, subject to dynamics
sphere((0, 0, 0.8), .5),
sphere((0, 0.8, 0), .5)
]
)
# FUNCTIONAL COMPONENTS
# simulation loop -- see presentation for the explanation
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
),
# Apply gravity force to particles. damping: numerical dissipation of
energy.
NewtonIntegrator(gravity=(0, 0, 0), damping=0.1)
]
# set timestep to a fraction of the critical timestep
# the fraction is very small, so that the simulation is not too fast
# and the motion can be observed
O.dt = .5e-4 * PWaveTimeStep()
O.step()
print O.forces.f(0)
a=[]
for i in O.bodies[0].intrs():
f = i.phys.normalForce + i.phys.shearForce
a.append(f)
resultF=Vector3(0,0,0)
for i in range(len(a)):
resultF=resultF+a[i]
print resultF
########### MWE2 which can lead to different results###########
readParamsFromTable(rMean=.05, rRelFuzz=.3, maxLoad=1e6, minLoad=1e4)
# make rMean, rRelFuzz, maxLoad accessible directly as variables later
from yade.params.table import *
# create box with free top, and ceate loose packing inside the box
from yade import pack, plot
O.bodies.append(geom.facetBox((.5, .5, .5), (.5, .5, .5), wallMask=31))
sp = pack.SpherePack()
sp.makeCloud((0, 0, 0), (1, 1, 1), rMean=rMean, rRelFuzz=rRelFuzz)
sp.toSimulation()
O.engines = [
ForceResetter(),
# sphere, facet, wall
InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb(),
Bo1_Wall_Aabb()]),
InteractionLoop(
# the loading plate is a wall, we need to handle sphere+sphere,
sphere+facet, sphere+wall
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom(),
Ig2_Wall_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
),
NewtonIntegrator(gravity=(0, 0, -9.81), damping=0.5),
# the label creates an automatic variable referring to this engine
# we use it below to change its attributes from the functions called
PyRunner(command='checkUnbalanced()', realPeriod=2, label='checker'),
]
O.dt = .5 * PWaveTimeStep()
# the following checkUnbalanced, unloadPlate and stopUnloading functions are
all called by the 'checker'
# (the last engine) one after another; this sequence defines progression of
different stages of the
# simulation, as each of the functions, when the condition is satisfied,
updates 'checker' to call
# the next function when it is run from within the simulation next time
# check whether the gravity deposition has already finished
# if so, add wall on the top of the packing and start the oedometric test
def checkUnbalanced():
# at the very start, unbalanced force can be low as there is only few
contacts, but it does not mean the packing is stable
if O.iter < 5000:
return
# the rest will be run only if unbalanced is < .1 (stabilized packing)
if unbalancedForce() > .1:
return
# add plate at the position on the top of the packing
# the maximum finds the z-coordinate of the top of the topmost particle
O.bodies.append(wall(max([b.state.pos[2] + b.shape.radius for b in
O.bodies if isinstance(b.shape, Sphere)]), axis=2, sense=-1))
global plate # without this line, the plate variable would only exist
inside this function
plate = O.bodies[-1] # the last particles is the plate
# Wall objects are "fixed" by default, i.e. not subject to forces
# prescribing a velocity will therefore make it move at constant
velocity (downwards)
plate.state.vel = (0, 0, -.1)
# start plotting the data now, it was not interesting before
O.engines = O.engines + [PyRunner(command='addPlotData()',
iterPeriod=200)]
# next time, do not call this function anymore, but the next one
(unloadPlate) instead
checker.command = 'unloadPlate()'
def unloadPlate():
# if the force on plate exceeds maximum load, start unloading
if abs(O.forces.f(plate.id)[2]) > maxLoad:
plate.state.vel *= -1
# next time, do not call this function anymore, but the next
one (stopUnloading) instead
checker.command = 'stopUnloading()'
def stopUnloading():
if abs(O.forces.f(plate.id)[2]) < minLoad:
# O.tags can be used to retrieve unique identifiers of the
simulation
# if running in batch, subsequent simulation would overwrite
each other's output files otherwise
# d (or description) is simulation description (composed of
parameter values)
# while the id is composed of time and process number
plot.saveDataTxt(O.tags['d.id'] + '.txt')
O.pause()
def addPlotData():
if not isinstance(O.bodies[-1].shape, Wall):
plot.addData()
return
Fz = O.forces.f(plate.id)[2]
plot.addData(Fz=Fz, w=plate.state.pos[2] - plate.state.refPos[2],
unbalanced=unbalancedForce(), i=O.iter)
O.run(10000,1)
waitIfBatch()
print O.forces.f(20)
a=[]
for i in O.bodies[20].intrs():
f = i.phys.normalForce + i.phys.shearForce
a.append(f)
resultF=Vector3(0,0,0)
for i in range(len(a)):
resultF=resultF+a[i]
print resultF
Thanks
Leonard
[1]https://yade-dev.gitlab.io/trunk/tutorial-examples.html
--
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 : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
