New question #706256 on Yade:
https://answers.launchpad.net/yade/+question/706256
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=200, 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)
# Generate sphere pack
O.materials.append(FrictMat(young=1E6, poisson=20000/1E6, density=10,
label='agg'))
sp = pack.SpherePack()
sp.makeCloud((-0.05,-0.05,0),
(0.05,0.05,0.09),rMean=0.01575,rRelFuzz=0,periodic=False,num=6)
sp2_1 = pack.SpherePack()
sp2_1.makeCloud((-0.03,-0.03,0.1),
(0.03,0.03,0.15),rMean=0.01175,rRelFuzz=0,periodic=False,num=7)
sp2_2 = pack.SpherePack()
sp2_2.makeCloud((-0.03,-0.03,0.16),
(0.03,0.03,0.2),rMean=0.011,rRelFuzz=0,periodic=False,num=8)
sp2_3 = pack.SpherePack()
sp2_3.makeCloud((-0.03,-0.03,0.21),
(0.03,0.03,0.25),rMean=0.01025,rRelFuzz=0,periodic=False,num=8)
sp3_1 = pack.SpherePack()
sp3_1.makeCloud((-0.03,-0.03,0.26),
(0.03,0.03,0.3),rMean=0.0083125,rRelFuzz=0,periodic=False,num=31)
sp3_2 = pack.SpherePack()
sp3_2.makeCloud((-0.03,-0.03,0.31),
(0.03,0.03,0.35),rMean=0.007125,rRelFuzz=0,periodic=False,num=31)
sp3_3 = pack.SpherePack()
sp3_3.makeCloud((-0.03,-0.03,0.36),
(0.03,0.03,0.4),rMean=0.0059375,rRelFuzz=0,periodic=False,num=32)
sp4_1 = pack.SpherePack()
sp4_1.makeCloud((-0.03,-0.03,0.41),
(0.03,0.03,0.45),rMean=0.0041525,rRelFuzz=0,periodic=False,num=25)
sp4_2 = pack.SpherePack()
sp4_2.makeCloud((-0.03,-0.03,0.46),
(0.03,0.03,0.5),rMean=0.003555,rRelFuzz=0,periodic=False,num=26)
sp4_3 = pack.SpherePack()
sp4_3.makeCloud((-0.03,-0.03,0.51),
(0.03,0.03,0.55),rMean=0.0029575,rRelFuzz=0,periodic=False,num=26)
sp5_1 = pack.SpherePack()
sp5_1.makeCloud((-0.03,-0.03,0.56),
(0.03,0.03,0.6),rMean=0.002065,rRelFuzz=0,periodic=False,num=25)
sp5_2 = pack.SpherePack()
sp5_2.makeCloud((-0.03,-0.03,0.61),
(0.03,0.03,0.65),rMean=0.00177,rRelFuzz=0,periodic=False,num=26)
sp5_3 = pack.SpherePack()
sp5_3.makeCloud((-0.03,-0.03,0.66),
(0.03,0.03,0.7),rMean=0.001475,rRelFuzz=0,periodic=False,num=26)
sp6_1 = pack.SpherePack()
sp6_1.makeCloud((-0.03,-0.03,0.56),
(0.03,0.03,0.6),rMean=0.001035,rRelFuzz=0,periodic=False,num=155)
sp6_2 = pack.SpherePack()
sp6_2.makeCloud((-0.03,-0.03,0.61),
(0.03,0.03,0.65),rMean=0.00089,rRelFuzz=0,periodic=False,num=155)
sp6_3 = pack.SpherePack()
sp6_3.makeCloud((-0.03,-0.03,0.66),
(0.03,0.03,0.7),rMean=0.000745,rRelFuzz=0,periodic=False,num=156)
sp7_1 = pack.SpherePack()
sp7_1.makeCloud((-0.03,-0.03,0.71),
(0.03,0.03,0.75),rMean=0.000525,rRelFuzz=0,periodic=False,num=155)
sp7_2 = pack.SpherePack()
sp7_2.makeCloud((-0.03,-0.03,0.76),
(0.03,0.03,0.8),rMean=0.00045,rRelFuzz=0,periodic=False,num=155)
sp7_3 = pack.SpherePack()
sp7_3.makeCloud((-0.03,-0.03,0.81),
(0.03,0.03,0.85),rMean=0.000375,rRelFuzz=0,periodic=False,num=156)
sp8_1 = pack.SpherePack()
sp8_1.makeCloud((-0.03,-0.03,0.86),
(0.03,0.03,0.9),rMean=0.0002625,rRelFuzz=0,periodic=False,num=155)
sp8_2 = pack.SpherePack()
sp8_2.makeCloud((-0.03,-0.03,0.91),
(0.03,0.03,0.95),rMean=0.000225,rRelFuzz=0,periodic=False,num=155)
sp8_3 = pack.SpherePack()
sp8_3.makeCloud((-0.03,-0.03,0.96),
(0.03,0.03,1),rMean=0.0001875,rRelFuzz=0,periodic=False,num=156)
sp9 = pack.SpherePack()
sp9.makeCloud((-0.03,-0.03,0.01),
(0.03,0.03,0.05),rMean=0.000075,rRelFuzz=0,periodic=False,num=8000)
sp10 = pack.SpherePack()
sp10.makeCloud((-0.03,-0.03,0.04),
(0.03,0.03,1),rMean=0.00005,rRelFuzz=0,periodic=False,num=80000)
# add the sphere pack to the simulation
sp.toSimulation(material='agg')
sp2_1.toSimulation(material='agg')
sp2_2.toSimulation(material='agg')
sp2_3.toSimulation(material='agg')
sp3_1.toSimulation(material='agg')
sp3_2.toSimulation(material='agg')
sp3_3.toSimulation(material='agg')
sp4_1.toSimulation(material='agg')
sp4_2.toSimulation(material='agg')
sp4_3.toSimulation(material='agg')
sp5_1.toSimulation(material='agg')
sp5_2.toSimulation(material='agg')
sp5_3.toSimulation(material='agg')
sp6_1.toSimulation(material='agg')
sp6_2.toSimulation(material='agg')
sp6_3.toSimulation(material='agg')
sp7_1.toSimulation(material='agg')
sp7_2.toSimulation(material='agg')
sp7_3.toSimulation(material='agg')
sp8_1.toSimulation(material='agg')
sp8_2.toSimulation(material='agg')
sp8_3.toSimulation(material='agg')
sp9.toSimulation(material='agg')
sp10.toSimulation(material='agg')
# 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.2),
NewtonIntegrator(gravity=(0,0,-1000), damping=0.4),
]
# Define simulation duration and time step
O.dt = 0.5e-6
O.run(10, wait=True)
for body in O.bodies:
if not isinstance(body.shape, Sphere):
continue
if body.shape.radius == 0.01575: #SP
body.shape.color = (0,0,0)
if body.shape.radius == 0.01175 :
body.shape.color = (0,0,1)
if body.shape.radius == 0.011 :
body.shape.color = (0,0,1)
if body.shape.radius == 0.01025 :
body.shape.color = (0,0,1)
if body.shape.radius == 0.00831250 : #SP3
body.shape.color = (0,1,0)
if body.shape.radius == 0.007125 :
body.shape.color = (0,1,0)
if body.shape.radius == 0.0059375 :
body.shape.color = (0,1,0)
if body.shape.radius == 0.0041525: #SP4
body.shape.color = (1,0,0)
if body.shape.radius == 0.003555 :
body.shape.color = (1,0,0)
if body.shape.radius == 0.0029575 :
body.shape.color = (1,0,0)
if body.shape.radius == 0.002065: #SP5
body.shape.color = (1,1,0)
if body.shape.radius == 0.00177 :
body.shape.color = (1,1,0)
if body.shape.radius == 0.001475 :
body.shape.color = (1,1,0)
if body.shape.radius == 0.001035: #SP6
body.shape.color = (1,1,1)
if body.shape.radius == 0.00089 :
body.shape.color = (1,1,1)
if body.shape.radius == 0.000745 :
body.shape.color = (1,1,1)
if body.shape.radius == 0.000525: #SP7
body.shape.color = (1,0.5,1)
if body.shape.radius == 0.00045 :
body.shape.color = (1,0.5,1)
if body.shape.radius == 0.000375 :
body.shape.color = (1,0.5,1)
if body.shape.radius == 0.0002625: #SP8
body.shape.color = (1,1,0.5)
if body.shape.radius == 0.000225 :
body.shape.color = (1,1,0.5)
if body.shape.radius == 0.0001875 :
body.shape.color = (1,1,0.5)
if body.shape.radius == 0.000075: #SP9
body.shape.color = (0.5,1,0.5)
if body.shape.radius == 0.00005: #SP10
body.shape.color = (0.5,0.5,1)
def checkUnbalanced():
if unbalancedForce() < .05:
O.pause()
plot.saveDataTxt('bbb.txt.bz2')
# plot.saveGnuplot('bbb') is also possible
# 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:
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, "%")
When I try to simulate I don't understand why the particles fall under the
container?
--
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
