yade-users team mailing list archive

[Huan <question706232@xxxxxxxxxxxxxxxxxxxxx>]

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

import random
import math
from yade import geom, pack, utils, ymport
from yade import export
import numpy as np

# Define cylinder parameters
center = (0, 0, 0)
radius = 0.102
height = 0.064

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

# define material properties
mat = 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 sphere parameters and percentages
radii = [0.01575, 0.01175, 0.011, 0.01025, 0.0083125, 0.007125, 0.0059375, 0.0041525, 0.003555, 0.0029575, 0.002065, 0.00177, 0.001475, 0.001035, 0.00089, 0.000745, 0.000525, 0.00045, 0.000375, 0.0002625, 0.000225, 0.0001875, 0.000075]
numSphere = [6, 7, 8, 8, 31, 31, 32, 25, 26, 26, 25, 26, 26, 155, 155, 156, 155, 156, 156, 155, 156, 156, 8000]

# Generate sphere pack
sp = pack.SpherePack()
for r, num in zip(radii, numSphere):
    sp.makeCloud((-0.055,-0.055,0.35), (0.055,0.055,0.01), rMean=r, rRelFuzz=0, num=num)

# add the sphere pack to the simulation
sp.toSimulation()

# Define gravity engine
O.engines = [    ForceResetter(),    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom6D(), Ig2_Facet_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_FrictPhys()],
        [Law2_ScGeom_FrictPhys_CundallStrack()]
    ),
    GlobalStiffnessTimeStepper(defaultDt=0.0001, timestepSafetyCoefficient=0.8),
    NewtonIntegrator(gravity=(0,0,-9.81), damping=0.4),
]

# Define simulation duration and time step
O.dt = 1e-5
O.run(1000, wait=True)

# Count the number of spheres in the cylinder
num_spheres_in_cylinder = 0
for body in O.bodies:
    if body.id == 0:  # skip the cylinder
        continue
    sphere_center = body.state.pos
    dist_to_axis = math.sqrt(sphere_center[0]**2 + sphere_center[1]**2)
    if dist_to_axis <= radius and sphere_center[2] >= center[2] and sphere_center[2] <= center[2] + height:
        num_spheres_in_cylinder += 1

# Calculate volume of spheres in cylinder
volume_of_spheres_in_cylinder = 0
for body in O.bodies:
    if body.id == 0:  # skip the cylinder
        continue
    sphere_center = body.state.pos
    dist_to_axis = math.sqrt(sphere_center[0]**2 + sphere_center[1]**2)
    if dist_to_axis <= radius and sphere_center[2] >= center[2] and sphere_center[2] <= center[2] + height:
        if isinstance(body.shape, Sphere):  # check if body is a sphere
            sphere_volume = (4/3) * math.pi * body.shape.radius**3
            volume_of_spheres_in_cylinder += sphere_volume

# Calculate volume of cylinder
volume_of_cylinder = (math.pi * radius**2) * height

# Calculate porosity
porosity = (volume_of_cylinder - volume_of_spheres_in_cylinder) / volume_of_cylinder
porosity_percent = porosity * 100

# Print results
print("Number of spheres in cylinder:", num_spheres_in_cylinder)
print("Volume of spheres in cylinder:", volume_of_spheres_in_cylinder)
print("Volume of cylinder:", volume_of_cylinder)
print("Porosity:", porosity)
print("Porosity:", porosity_percent, "%")

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