Is it possible to add hue to each radii sphere radius? Also, after simulation is my calculation wrong?

Asked by Huan

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

# Define cylinder parameters
center = (0, 0, 0)
radius = 0.08
height = 0.10

# create cylindrical body with radius 0.08 m and height 0.0605 m
cylinder = yade.geom.facetCylinder(center=center, radius=radius, height=height, segmentsNumber=80, wallMask=6)
O.bodies.append(cylinder)

# Define sphere parameters and percentages
radii = [0.000075, 0.0006, 0.00236, 0.00475, 0.0095, 0.0125]
percentages = [0.04, 0.06, 0.05, 0.35, 0.2, 0.3]

# Generate sphere pack
sp = pack.SpherePack()
for r, ratio in zip(radii, percentages):
   num_spheres = int(ratio * 1000)
   sp.makeCloud((-0.055,-0.055,0.35), (0.055,0.055,0.01), rMean=r, rRelFuzz=0, num=num_spheres)

# 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)

# 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, "%")

I have two question
1. Is it possible to add color to each radii sphere radius?
2. After simulation, the sphere seems to pack well in the container, but my porosity after calculation if above 80+%
did I have mistake in my calculation process?

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
Last query:
Last reply:
Revision history for this message
Best Jan Stránský (honzik) said :
#1

Hello,

> I have two question

next time, please open separate questions for separate problems ([1], point 5)

> 1. Is it possible to add color to each radii sphere radius?

yes. e.g.:
###
...
sp.toSimulation()
rMax = max(b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere))
rMin = min(b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere))
def radius2color(r,rMin,rMax):
    factor = (r-rMin)/(rMax-rMin)
    r = factor
    g = 0
    b = 1 - factor
    return (r,g,b)
for b in O.bodies:
    if not isinstance(b.shape,Sphere):
        continue
    b.shape.color = radius2color(b.shape.radius,rMin,rMax)
###

> O.run(1000)
> 2. After simulation, the sphere seems to pack well in the container, but my porosity after calculation if above 80+% did I have mistake in my calculation process?

The problem is that what you print is at the very beginning of the simulation, not "after calculation"
As mentioned in previous question, DO NOT use O.run(N) if you have code after it. Then the code is executed just after the O.run command, i.e. during the first time steps of the simulation.

Use
###
O.run(N,wait=True)
###
or
###
O.run(N)
O.wait()
###

If I use
O.run(20000,wait=True)
instead, I get porosity below 80%, although still very high.
Maybe some boundary effects..

Cheers
Jan

[1] https://www.yade-dem.org/wiki/Howtoask

Revision history for this message
Huan (huan-liu) said (last edit ):
#2

Thanks Jan Stránský, that solved my question.
Sorry I didn't know about the rule, but I need to modify some of my coding I think porosity isn't properly calculated