# Force applied from disk to particle

Hi,

I compacting spheres in cylinder, have wall as upper punch and a disk as lower punch. How do I calculate the force from disk applied on the particles.

For the walls I calculate the applied force as abs(O.forces.f(UP.id)[2])

For the disk I tried using:
for n in lower_punch:
body= O.bodies[n]
force_lp=abs(O.forces.f(body.id)[2])

That gives me a small number, which seems to be wrong.

Best regards,
Mithu

Here is my code:

from __future__ import print_function
from yade import utils, plot, timing
import pandas as pd
import numpy as np
from PIL import Image
from scipy.interpolate import interp1d

o = Omega()

# Physical parameters
fr = 0.41
rho = 1561
Diameter = 7.9e-5
D=Diameter
r1 = Diameter/2
#r2 = Diameter/2
k1 = 400
kp = 10.0*k1
kc = k1 * 0.1
ks = k1 * 0.1
DeltaPMax = Diameter/3.0
Chi1 = 0.34

o.dt = 1.0e-7
particleMass = 4.0/3.0*math.pi*r1*r1*r1*rho
m_tot=4.5e-06
no_p=m_tot/particleMass

Vi1 = math.sqrt(k1/particleMass)*DeltaPMax*Chi1
PhiF1=0.999
#PhiF1 = DeltaPMax*(kp-k1)*(r1+r2)/(kp*2*r1*r2)

Cyl_height=0.007

Comp_press_up= 0.6e8
Comp_force_up=Comp_press_up*cross_area

Comp_press_lp= 0.6e8
Comp_force_lp=Comp_press_lp*cross_area

compression_data_save=[]
sc_por=0.1
#*************************************

mat1 = O.materials.append(LudingMat(frictionAngle=fr, density=rho, k1=k1, kp=kp, ks=ks, kc=kc, PhiF=PhiF1, G0 = 0.0))

# Spheres for compression and walls
sp=pack.SpherePack()
sp.makeCloud((-7.5*Diameter,-7.5*Diameter,-44*Diameter),(7.5*Diameter,7.5*Diameter,44.0*Diameter), rMean=Diameter/2.0,rRelFuzz=0.18,num=round(no_p))
sp.toSimulation(material=mat1)

o.engines = [
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.05),
Bo1_Wall_Aabb(),
Bo1_Facet_Aabb()
]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.05),
Ig2_Facet_Sphere_ScGeom(),
Ig2_Wall_Sphere_ScGeom()],
[Ip2_LudingMat_LudingMat_LudingPhys()],
[Law2_ScGeom_LudingPhys_Basic()]
),
NewtonIntegrator(damping=0.1, gravity=[0, 0, -9.81]),
PyRunner(command='checkForce()', realPeriod=1, label="fCheck"),
#DeformControl(label="DefControl")
]

def checkForce():
if O.iter < 1000:
return
#if unbalancedForce() > 0.8:
#return
O.bodies.append(wall((-Cyl_height/2)+utils.aabbDim()[2]+3*Diameter, axis=2, sense=-1, material=mat1))
global UP
UP = O.bodies[-1]
UP.state.vel = (0, 0, -5)
global lower_punch
for n in lower_punch:
body= O.bodies[n]
body.state.vel = (0,0,5)

for n in lower_punch:
body= O.bodies[n]
force_lp=abs(O.forces.f(body.id)[2])
if ((abs(O.forces.f(UP.id)[2]) > Comp_force_up) and (force_lp > Comp_force_lp)):
UP.state.vel *= -1
#LP.state.vel *= -1

for n in lower_punch:
body= O.bodies[n]
force_lp=abs(O.forces.f(body.id)[2])
if ((abs(O.forces.f(UP.id)[2])==0) and (force_lp==0)):
o.pause()
compression_pressure=((abs(O.forces.f(UP.id)[2])+force_lp)/(cross_area))*1e-6
porosity=voxelPorosity(resolution=200,start=(utils.aabbExtrema()[0]+(sc_por*D,sc_por*D,sc_por*D)),end=(utils.aabbExtrema()[1]-(sc_por*D,sc_por*D,sc_por*D)))
data_to_save=[porosity,compression_pressure]
compression_data_save.append(data_to_save)

## Question information

Language:
English Edit question
Status:
Solved
For:
Assignee:
No assignee Edit question
Solved by:
Mithushan Soundaranathan
Solved:
Last query:
 Revision history for this message Karol Brzezinski (kbrzezinski) said on 2021-09-01: #1

Hi!

It seems that you have forgotten about summing the forces from separate facets ;)
Try this:
##############
force_lp = 0
for n in lower_punch:
body = O.bodies[n]
force_lp = force_lp + abs(O.forces.f(body.id)[2])
##############

Best wishes,
Karol

 Revision history for this message Jan Stránský (honzik) said on 2021-09-02: #2

or, in a "more pythonic" way:

force_lp = sum(abs(O.forces.f(n)[2]) for n in lower_punch)

Also, is the abs necessary? Should not..

Cheers
Jan

 Revision history for this message Mithushan Soundaranathan (mithushan93) said on 2021-09-11: #3

Hi Both,

Thank you very much, it was very helpful for me.

Best,
Mithu