yade-users team mailing list archive

[Mithushan Soundaranathan <question697108@xxxxxxxxxxxxxxxxxxxxx>]

New question #697108 on Yade:
https://answers.launchpad.net/yade/+question/697108

Hi,
I have compacted spherical particles in a cylinder, and the particle are overlapping.  I am struggling to measure the porosity of my compacys. I have tried utils.porosity, voxelporosity and also tried using the usual porosity Eq. (1- (rho_compacts/rho_t)). They are giving me either negative values upto (-60%) or to high values (without accounting for the overlap). I have looked through the forum for similar problems, I have not any solution which applicable for my simulation. I was wondering if anyone have experienced similar issues and how they solved it.

Best,
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


o = Omega()

# Physical parameters
fr        = 0.54
rho       = 1050
Diameter  = 0.0012
D=Diameter
r1        = Diameter/2
#r2        = Diameter/2
k1        = 126000
kp        = 12.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

Vi1 = math.sqrt(k1/particleMass)*DeltaPMax*Chi1
PhiF1=0.999


Tab_rad=0.005
Cyl_height=0.02

Comp_press1=1.4e8
Comp_force1=Comp_press1*cross_area

# 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
sp=pack.SpherePack()
sp.makeCloud((-3.0*Diameter,-3.0*Diameter,-8*Diameter),(3.0*Diameter,3.0*Diameter,28.0*Diameter), rMean=Diameter/2.0,num=527)
#cyl = pack.inCylinder((0,0,0),(0,0,40*Diameter),40*Diameter-0.006)
#sp = pack.filterSpherePack(cyl,sp,True,material=mat1)
sp.toSimulation(material=mat1)



######################################################################
#O.bodies.append(geom.facetBox((0,0,0), (4.0*Diameter,4.0*Diameter,4.0*Diameter), wallMask=63-32, 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]),
  #VTKRecorder(fileName='vtk-',recorders=['all'],iterPeriod=10000),
  PyRunner(command='checkForce()', realPeriod=1, label="fCheck"),
  #DeformControl(label="DefControl")
]

def checkForce():
    if O.iter < 1000000:
        return
    timing.reset()
    if unbalancedForce() > 0.2:
        return
    # add plate at upper box side

    highSphere = 0.0
    for b in O.bodies:
        if highSphere < b.state.pos[2] and isinstance(b.shape, Sphere):
            highSphere = b.state.pos[2]
        else:
            pass

    O.bodies.append(wall(highSphere+0.5*Diameter, axis=2, sense=-1, material=mat1))
    global plate
    plate = O.bodies[-1]  
    plate.state.vel = (0, 0, -.3)
    fCheck.command = 'unloadPlate()'


def unloadPlate():
    if abs(O.forces.f(plate.id)[2]) > Comp_force1:
        plate.state.vel *= -1
        fCheck.command = 'stopUnloading()'


def stopUnloading():
    if abs(O.forces.f(plate.id)[2]) == 0:
       O.pause()



    
    





    





-- 
You received this question notification because your team yade-users is
an answer contact for Yade.