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.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
