Obtaining Volume Fraction of Filled Cylinder

Asked by Kevin on 2020-04-26

Hey guys! So I have been having trouble finding a way in which to obtain the volume fraction for a filled coaxial cylinder model of spheres using utils.getDepthProfiles. I imported a quarter cylinder model via solidworks (to minimize computational complexity) and filled the outer coaxial cylinder with thousands of spheres with varying radii. The objective here is to use the utils.getDepthProfile for a portion of the cylinder volume and I am unsure how to print out the resultant volume fraction. YadeBook helps me with the parameters (the volume of the right cylinder is the first float parameter that I had calculated previously), but I am unsure about how to use the latter two parameters, or better yet, how to use the function to extract the volume fraction at all (I am a novice at programming). Under "Tracking Volume Fraction" is where I am having trouble implementing or finding a spot to extract the overall volume fraction after filling the model. How should I approach this based on this code sample?

Thank you for any help or advice given!

#Module Setup
from yade import ymport, plot
import time

#Constants
utils.readParamsFromTable(descriptionIn = 'fill',
  frIn = 0.5, enIn=0.01, etIn=0.01, tcIn=0.0001, # Frictionangle (rad), strength-parameter
  rhoIn = 2300.0, # Density
  dumpVTKIn = 100000 # Periods of dumps in iterations
)

#Sweep Parameters using Table
from yade.params.table import
import shutil

try:
  shutil.rmtree('step1')
except OSError:
  pass
os.mkdir('step1')

folderNameBase = 'step1/' + str(descriptionIn)
folderName = folderNameBase
os.mkdir(folderNameBase)

#Initialize
o = Omega()
o.dt = 0.05*1E-7

#Define Material Properties of Spheres
mat1 = o.materials.append(ViscElMat(frictionAngle=frIn, density=rhoIn, en=enIn, et=etIn))

#Build Geometry
cylinder_outer=o.bodies.append(geom.facetCylinder(center=(0.,0.,0.),radius=3.5e-3,height=5e-3,
    radiusTopInner=1.5e-3,radiusBottomInner=1.5e-3,segmentsNumber=80,wallMask=6,closeGap=True))
cylinder_inner=o.bodies.append(geom.facetCylinder(center=(0.,0.,0.),radius=1.5e-3,height=5e-3,
    segmentsNumber=80,wallMask=7,closeGap=True))
cylinder=o.bodies.append(ymport.stl('coax_2mm_coneb.SLDPRT'))

o.engines = [
    #Following Default Workflow
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()],sortThenCollide=False,verletDist=-0.5,label='collider'),
# SpatialQuickSortCollider(),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_FrictPhys()],
        [Law2_ScGeom_FrictPhys_CundallStrack()],
    ),
    NewtonIntegrator(damping=0.5,gravity=Vector3(0,-90.81,0), label='newt'),
    #Particles will drop from a circular plane
    CircularFactory(maxParticles=5000, #Number of generated particles
        center=(0.0,2.5e-3,0.0),
        radius = (3.5e-3),
        PSDcalculateMass=False,
        PSDsizes = [125e-6,
                    146e-6,
                    167e-6,
                    188e-6,
                    209e-6,
                    230e-6,
                    251e-6,
                    272e-6,
                    293e-6,
                    314e-6,
                    335e-6,
                    356e-6,
                    376e-6,
                    397e-6,
                    418e-6,
                    439e-6,
                    460e-6],
        PSDcum = [0.0250,
                    0.0328,
                    0.1197,
                    0.3573,
                    0.6326,
                    0.7687,
                    0.8431,
                    0.8929,
                    0.9206,
                    0.9420,
                    0.9609,
                    0.9666,
                    0.9748,
                    0.9824,
                    0.9868,
                    0.9905,
                    0.9943,
                    ],
        vMin=1e-1,vMax=1e-0,vAngle=0,massFlowRate=10e-4,
        normal=(0.0,-1.0,0.0),label='factorySpheres',materialId=mat1,
        stopIfFailed=False,silent=True
    ),
    #Delete Particles that fly out of domain
    DomainLimiter(lo=(-4e-3,-1e-3,-4e-3),hi=(4e-3,4e-3,4e-3),iterPeriod=5000),
    #Shake the domain to help settling
    HarmonicMotionEngine(A=[0,0,0.1e-4],f=[0,0,5000],ids=cylinder),
    #Record output to VTK Files, Spheres,Interactions,Facets
    VTKRecorder(iterPeriod=dumpVTKIn,fileName=folderName+'/spheres-',
        recorders=['spheres','velocity','colors','intr','ids','mask','materialId','stress'], label='VTK_sphere'),
    VTKRecorder(iterPeriod=dumpVTKIn,fileName=folderName+'/facetB-',
        recorders=['facets','colors'],label='VTK_OD'),
    #FIXME Snapshot Taker
# SnapshotEngine(iterPeriod=dumpVTKIn,fileBase=folderName+'/spheres-',viewNo=0,label='snapshooter'),
    #Execute Python Commands
    PyRunner(iterPeriod=10000,initRun=True,command='stop_condition()'),
# PyRunner(iterPeriod=1000,command='addPlotData()')
    PyRunner(iterPeriod=5000000,command='o.save(\'save1.xml\')')
]

#Enable Track Energy - Can be used to determine if sim is done
o.trackEnergy=True

#Volume Fraction Tracking
from yade import utils
utils.getDepthProfiles(1.923e-7,40,1e-4,0)

#Run Yade Simulation
o.stopAtIter=5000001 #Maximum Iterations
o.run()

def stop_condition():
    o.dt=0.5*utils.PWaveTimeStep() #auto timesetup
    print('iter,time,dt,iter/sec=',o.iter,o.iter*o.dt,o.dt,o.speed) #print to screen iteration

Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
2020-04-26
Last reply:
2020-04-30
Jérôme Duriez (jduriez) said : #1

I'm not really familiar with that function and am even less sure I understood your question, but :

if getDepthProfiles() returns a tuple of "local" quantities, the overall quantity could probably be obtained averaging the local ones.

Did you try averaging the elements of the returned tuple (with appropriate weights ?)

PS : if the question is about getDepthProfiles(), why not include it in the title ? I could help reduce confusion..

Can you help with this problem?

Provide an answer of your own, or ask Kevin for more information if necessary.

To post a message you must log in.