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
from yade import pack
import pandas as pd
import numpy as np
from PIL import Image
from yade import pack, export
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)

Tab_rad=0.001
Cyl_height=0.007
cross_area=math.pi*(Tab_rad**2)

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
#*************************************

# Add material
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)
walls=O.bodies.append(yade.geom.facetCylinder((0,0,0),radius=Tab_rad,height=Cyl_height,segmentsNumber=20,wallMask=6,material=mat1))

# Add engines
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
    lower_punch= O.bodies.append(geom.facetCylinder((0,0,(-Cyl_height/2)),Tab_rad-.00001,0,segmentsNumber=50,wallMask=1))
    for n in lower_punch:
        body= O.bodies[n]
        body.state.vel = (0,0,5)
    fCheck.command = 'unloadPlate()'

def unloadPlate():
    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
        fCheck.command = 'stopUnloading()'

def stopUnloading():
    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)