# How to get the resultant force on a particle

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)

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)
# 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.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###########

# make rMean, rRelFuzz, maxLoad accessible directly as variables later

# create box with free top, and ceate loose packing inside the box
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 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 + 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

plate.state.vel *= -1

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

if not isinstance(O.bodies[-1].shape, Wall):
return
Fz = O.forces.f(plate.id)
plot.addData(Fz=Fz, w=plate.state.pos - plate.state.refPos, unbalanced=unbalancedForce(), i=O.iter)

O.run(10000,1)

waitIfBatch()

print O.forces.f(20)

a=[]
for i in O.bodies.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

## Question information

Language:
English Edit question
Status:
Solved
For:
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
Last query:
 Revision history for this message  Jan Stránský (honzik) said on 2021-11-16: #1

Hello,

there are two points:
1) NewtonInegrator.gravity. It is not included in O.forces.f
2) i.phys.normalForce and i.phys.shearForce are numbers/vectors SAME for both interacting bodies, but clearly they must be taken into account oppositely for each body (Newton's 3rd law of action and reaction)

Below is corrected code. For the second example, replace 0 with 20 (body id)
###
NewtonIntegrator(...,label="integrator")
...

print O.forces.f(0) + newton.gravity * O.bodies.state.mass

a=[]
for i in O.bodies.intrs():
f = i.phys.normalForce + i.phys.shearForce
if O.bodies.id == i.id1: f = -f # HERE
a.append(f)

resultF=Vector3(0,0,0)
for i in range(len(a)):
resultF=resultF+a[i]
resultF += integrator.gravity * O.bodies.state.mass # HERE
###

Cheers
Jan

PS: Not important notes:

---
> for i in range(len(a)):
> resultF=resultF+a[i]

for f in a:
resultF=resultF+f

or even more "pythonic" without (explicit) loop:
resultF = sum(a, Vector3.Zero)

---
print O.forces.f(20)

a=[]
for i in O.bodies.intrs():

20 is repeated code/value, which is generally considered as bad practice (to change it, you need to be careful to change all occurrences, even more emphasized with the corrected code).
Better is to use a variable for this:
###
id = 20
print O.forces.f(id)
a=[]
for i in O.bodies[id].intrs():

 Revision history for this message Leonard (z2521899293) said on 2021-11-17: #2

Thanks Jan Stránský, that solved my question.