yade-users team mailing list archive

[Huan <question706332@xxxxxxxxxxxxxxxxxxxxx>]

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

import random
import math
from yade import geom, pack, utils
import pandas as pd

# Define cylinder parameters
center = (0, 0, 0)
cyl_radius = 0.102
cyl_height = 0.064

# create cylindrical body with radius 0.102 m and height 0.064 m
cylinder = geom.facetCylinder(center=center, radius=cyl_radius, height=cyl_height, segmentsNumber=80, wallMask=6)

# define material properties
mat = utils.PolyhedraMat()
mat.density = 2600  # kg/m^3
mat.young = 1E6  # Pa
mat.poisson = 20000 / 1E6
mat.frictionAngle = 0.6  # rad

# assign material to each body in the cylinder
for body in cylinder:
    body.bodyMat = mat

# add cylinder to simulation
O.bodies.append(cylinder)

# Define box dimensions
box_dimensions = (0.1, 0.1, 0.1)

# Create box at origin
box = geom.facetBox(center=(0,0,0), extents=box_dimensions, wallMask=31)

# Define material properties for box
mat = O.materials.append(FrictMat(density=2600, young=1E6, poisson=0.3, frictionAngle=0.6))

# Assign material to each facet of box
for body in box:
    body.bodyMat = mat

# Add box to simulation
O.bodies.append(box)

# define sphere parameters and number of spheres
radius = 0.01575
num = 6
radius1 = 0.01175
num1 = 7
radius2 = 0.011
num2 = 8
radius3 = 0.01025
num3 = 8
radius4 = 0.0083125
num4 = 31
radius5 = 0.007125
num5 = 31
radius6 = 0.0059375
num6 = 32
radius7 = 0.0041525
num7 = 25
radius8 = 0.003555
num8 = 26
radius9 = 0.0029575
num9 = 26
radius10 = 0.002065
num10 = 25
radius11 = 0.00177
num11 = 26
radius12 = 0.001475
num12 = 26
radius13 = 0.001035
num13 = 155
radius14 = 0.00089
num14 = 155
radius15 = 0.000745
num15 = 156
radius16 = 0.00045
num16 = 9735

## create empty sphere packing
sp = pack.SpherePack()
sp1 = pack.SpherePack()
sp2 = pack.SpherePack()
sp3 = pack.SpherePack()
sp4 = pack.SpherePack()
sp5 = pack.SpherePack()
sp6 = pack.SpherePack()
sp7 = pack.SpherePack()
sp8 = pack.SpherePack()
sp9 = pack.SpherePack()
sp10 = pack.SpherePack()
sp11 = pack.SpherePack()
sp12 = pack.SpherePack()
sp13 = pack.SpherePack()
sp14 = pack.SpherePack()
sp15 = pack.SpherePack()
sp16 = pack.SpherePack()

# generate randomly sphere
sp.makeCloud((-0.055,-0.055,0.0105), (0.055,0.055,0.0105), rMean=radius, rRelFuzz=0, num=num)
sp1.makeCloud((-0.055,-0.055,0.0380), (0.055,0.055,0.0380), rMean=radius1, rRelFuzz=0, num=num1)
sp2.makeCloud((-0.055,-0.055,0.0608), (0.055,0.055,0.0608), rMean=radius2, rRelFuzz=0, num=num2)
sp3.makeCloud((-0.055,-0.055,0.0820), (0.055,0.055,0.0820), rMean=radius3, rRelFuzz=0, num=num3)
sp4.makeCloud((-0.0605,-0.0605,0.1006), (0.0605,0.0605,0.1006), rMean=radius4, rRelFuzz=0, num=num4)
sp5.makeCloud((-0.0615,-0.0615,0.1160), (0.0615,0.0615,0.1160), rMean=radius5, rRelFuzz=0, num=num5)
sp6.makeCloud((-0.0605,-0.0605,0.1291), (0.0605,0.0605,0.1291), rMean=radius6, rRelFuzz=0, num=num6)
sp7.makeCloud((-0.0605,-0.0605,0.1392), (0.0605,0.0605,0.1392), rMean=radius7, rRelFuzz=0, num=num7)
sp8.makeCloud((-0.0605,-0.0605,0.1469), (0.0605,0.0605,0.1469), rMean=radius8, rRelFuzz=0, num=num8)
sp9.makeCloud((-0.0605,-0.0605,0.1534), (0.0605,0.0605,0.1534), rMean=radius9, rRelFuzz=0, num=num9)
sp10.makeCloud((-0.0605,-0.0605,0.1584), (0.0605,0.0605,0.1584), rMean=radius10, rRelFuzz=0, num=num10)
sp11.makeCloud((-0.0605,-0.0605,0.1622), (0.0605,0.0605,0.1622), rMean=radius11, rRelFuzz=0, num=num11)
sp12.makeCloud((-0.0605,-0.0605,0.1655), (0.0605,0.0605,0.1655), rMean=radius12, rRelFuzz=0, num=num12)
sp13.makeCloud((-0.0605,-0.0605,0.1680), (0.0605,0.0605,0.1680), rMean=radius13, rRelFuzz=0, num=num13)
sp14.makeCloud((-0.0605,-0.0605,0.1699), (0.0605,0.0605,0.1699), rMean=radius14, rRelFuzz=0, num=num14)
sp15.makeCloud((-0.0605,-0.0605,0.1715), (0.0605,0.0605,0.1715), rMean=radius15, rRelFuzz=0, num=num15)
sp16.makeCloud((-0.0605,-0.0605,0.1727), (0.0605,0.0605,0.1727), rMean=radius16, rRelFuzz=0, num=num16)

# add interaction of gravity of sphere and cylinder
O.engines = [ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]),
        InteractionLoop(
                # handle sphere+sphere and facet+sphere collisions
                [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()],
                [Ip2_FrictMat_FrictMat_FrictPhys()],
                [Law2_ScGeom_FrictPhys_CundallStrack()]
        ),
        NewtonIntegrator(gravity=(0,0,-9.81), damping=0.4),
]
O.dt = 0.1 * PWaveTimeStep()

# run the simulation for 1000 steps
O.run(1000, wait=True)

for body in O.bodies:
   if not isinstance(body.shape, Sphere): 
       continue
   if body.shape.radius == radius: #SP
       body.shape.color = (0,0,1) #blue
   if body.shape.radius == radius1: #SP1
       body.shape.color = (1,0,0) #red
   if body.shape.radius == radius2: #SP2
       body.shape.color = (0,1,0) #green
   if body.shape.radius == radius3: #SP3
       body.shape.color = (1,1,1) #white
   if body.shape.radius == radius4: #SP4
       body.shape.color = (1,1,0) #yellow
   if body.shape.radius == radius5: #SP5
       body.shape.color = (1,0.5,0) #orange
   if body.shape.radius == radius6: #SP6
       body.shape.color = (0.5,0,1) #purple
   if body.shape.radius == radius7: #SP7
       body.shape.color = (0.5,0.5,0.5) #gray
   if body.shape.radius == radius8: #SP8
       body.shape.color = (1,0.5,0.5) #pink
   if body.shape.radius == radius9: #SP9
       body.shape.color = (0,0.85,0.8) #turquoise
   if body.shape.radius == radius10: #SP10
       body.shape.color = (0.6,0.4,0.2) #brown
   if body.shape.radius == radius11: #SP11
       body.shape.color = (0.5,0,0) #maroon
   if body.shape.radius == radius12: #SP12
       body.shape.color = (0.5,0.5,0) #olive
   if body.shape.radius == radius13: #SP13
       body.shape.color = (1,0,1) #magenta
   if body.shape.radius == radius14: #SP14
       body.shape.color = (0,1,1) #cyan
   if body.shape.radius == radius15: #SP15
       body.shape.color = (1,0.84,0) #gold
   if body.shape.radius == radius16: #SP16
       body.shape.color = (0,0,0) #black


# calculate the volume of the packing
volume_packing = 0
num_spheres = 0
for b in O.bodies:
    if isinstance(b.shape, yade.wrapper.Sphere):
        volume_packing += 4/3 * math.pi * b.shape.radius**3
        num_spheres += 1

# calculate the volume of the cylinder
volume_cylinder = math.pi * cyl_radius**2 * cyl_height

# calculate the porosity and porosity percentage
porosity = (volume_cylinder - volume_packing) / volume_cylinder
porosity_percent = porosity * 100

print("Number of spheres:", "{:d}".format(num_spheres))
print("V Packing:", "{:e}".format(volume_packing))
print("V Cylinder:", "{:e}".format(volume_cylinder))
print("Porosity:", "{:.2f}".format(porosity))
print("Porosity:", "{:.2f}%".format(porosity_percent))

# print the number of spheres in each layer
print("Number of spheres in each layer:")
print(f"Layer 1: {len(sp)}")
print(f"Layer 2: {len(sp1)}")
print(f"Layer 3: {len(sp2)}")
print(f"Layer 4: {len(sp3)}")
print(f"Layer 5: {len(sp4)}")
print(f"Layer 6: {len(sp5)}")
print(f"Layer 7: {len(sp6)}")
print(f"Layer 8: {len(sp7)}")
print(f"Layer 9: {len(sp8)}")
print(f"Layer 10: {len(sp9)}")
print(f"Layer 11: {len(sp10)}")
print(f"Layer 12: {len(sp11)}")
print(f"Layer 13: {len(sp12)}")
print(f"Layer 14: {len(sp13)}")
print(f"Layer 15: {len(sp14)}")
print(f"Layer 16: {len(sp15)}")
print(f"Layer 17: {len(sp16)}")

# create a dictionary to store simulation data
data = {'Number of spheres': [num_spheres], 'V Packing': [volume_packing], 'V Cylinder': [volume_cylinder],
        'Porosity': [porosity], 'Porosity %': [porosity_percent]}

# create a DataFrame from the data
df = pd.DataFrame(data)

# append the data to an existing file or create a new file if it doesn't exist
with open('simulation_data.csv', 'a') as file:
    df.to_csv(file, header=file.tell() == 0, index=False)

1. Is it possible to only allows the particles to fall into the container without bouncing out after contact whether the contact might be against other spheres or the container?

In the simulation, I tried to generate a box and try to contain all the sphere into the cylinder, but it didn't work as I expected it to 

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