yade-users team mailing list archive

[Mithushan Soundaranathan <question696689@xxxxxxxxxxxxxxxxxxxxx>]

Question #696689 on Yade changed:
https://answers.launchpad.net/yade/+question/696689

Description changed to:
Hi,

Has anyone experience in measuring the porosity of particle compacts in
cylinder? I tried voxel porosity and I can only define domain size in
boxes, it is affected by the boundary and "air" on the compact edge.
Also is it possible define the domain size in cylinder for calculation
of voxel porosity?

Here is my code:
#!/usr/bin/env python
#encoding: ascii

# Testing of the Deformation Enginge with Luding Contact Law
# Modified Oedometric Test
# The reference paper [Haustein2017]


from __future__ import print_function
from yade import utils, plot, timing
from yade import pack
import pandas as pd

o = Omega()

# Physical parameters
fr        = 0.54
rho       = 1050
Diameter  = 0.0079
D=Diameter
r1        = Diameter
r2        = Diameter
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-5

particleMass = 4.0/3.0*math.pi*r1*r1*r1*rho

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

#surf_a=4*math.pi*(7.9e-5*7.9e-5)
#comp_stress=300e3*surf_a


#*************************************

# 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((-11.0*Diameter,-11.0*Diameter,-25*Diameter),(11*Diameter,11*Diameter,30.0*Diameter), rMean=Diameter/2.0,rRelFuzz=.18)
cyl = pack.inCylinder((0,0,0),(0,0,30*Diameter),11*Diameter-0.006)
sp = pack.filterSpherePack(cyl,sp,True,material=mat1)
sp.toSimulation(material=mat1)

##No of particles
count = 0
for b in O.bodies:
    if isinstance(b.shape, Sphere):
        count +=1

######################################################################
walls=O.bodies.append(yade.geom.facetCylinder((0,0,0),radius=11*Diameter,height=14*Diameter,segmentsNumber=20,wallMask=6,material=mat1))


df = pd.DataFrame(columns=['Cut_xy','Cut_z','Volume','Porosity'])
df.to_csv('PH101_rev.csv')
from csv import writer
def append_list_as_row(file_name, list_of_elem):
    # Open file in append mode
    with open(file_name, 'a+', newline='') as write_obj:
        # Create a writer object from csv module
        csv_writer = writer(write_obj)
        # Add contents of list as last row in the csv file
        csv_writer.writerow(list_of_elem)
# 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():
    # 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 < 20000:
        return
    # the rest will be run only if unbalanced is < .1 (stabilized packing)
    timing.reset()
    if unbalancedForce() > 0.4:
        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))
    # without this line, the plate variable would only exist inside this
    # function
    global plate
    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, -0.1)
    # start plotting the data now, it was not interesting before
    O.engines = O.engines + [PyRunner(command='addPlotData()', iterPeriod=1000)]
    # next time, do not call this function anymore, but the next one
    # (unloadPlate) instead
    fCheck.command = 'unloadPlate()'


def unloadPlate():
    # if the force on plate exceeds maximum load, start unloading
    # if abs(O.forces.f(plate.id)[2]) > 5e-2:
    if abs(O.forces.f(plate.id)[2]) > 83e3:
        plate.state.vel *= -1
        # next time, do not call this function anymore, but the next one
        # (stopUnloading) instead
        fCheck.command = 'stopUnloading()'


def stopUnloading():
    if abs(O.forces.f(plate.id)[2]) == 0:
        # 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')
        plot.saveDataTxt('data'+ O.tags['id'] +'.txt')
        #print(timing.stats())
        O.pause()
        Diameter_cut=[0.2,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5]
        height_tab=utils.aabbDim()[2]/(2*len(Diameter_cut))
        D=Diameter
        for i in range (0,len(Diameter_cut)):
            Cut_xy=Diameter_cut[i]
            Start_point=utils.aabbExtrema()[0]+(Cut_xy*D,Cut_xy*D,i*height_tab)
            End_point=utils.aabbExtrema()[1]-(Cut_xy*D,Cut_xy*D,i*height_tab)
            vol=(End_point[0]-Start_point[0])*(End_point[1]-Start_point[1])*(End_point[2]-Start_point[2])
            Porosity_rev= [voxelPorosity(resolution=1600,start=Start_point,end=End_point)]
            row_contents = [i,Start_point,End_point,vol,Porosity_rev]
            append_list_as_row('PH101_rev.csv', row_contents)
        
            


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


def defVisualizer():
   with open("data.dat","a") as f:
       for b in O.bodies:
           if isinstance(b.shape, Sphere):
               rData = "{x},{y},{z},{r},{w}\t".format(x = b.state.pos[0],
                                                y = b.state.pos[1],
                                                z = b.state.pos[2],
                                                r = b.shape.radius + b.state.dR,
                                                w = plate.state.pos[2]
                                               )
               f.write(rData)
       f.write("\n")


O.timingEnabled=True
O.run(1, True)
plot.plots={'w':('Fz', None)}
plot.plot()


Best,
Mithushan

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