Troble with porosity calculation
Hello,
I was able to create a sphere packing with gravity deposition and oedometric test. So, I create a calculation to calculate the porosity, but the result that I get is negative which doesn't make sense. The result is around 25.00%
The following is my code:
import random
import math
from yade import geom, pack, utils, plot, ymport
import pandas as pd
# Define material properties
youngModulus = 1e7
poissonRatio = 0.25
density = 2000
# Create material
material = O.materials.
# Define cylinder with funnel parameters
center = (0, 0, 0)
diameter = 0.102
height = 0.18
# create cylindrical body with radius 0.102 m and height 0.064 m
cylinder = geom.facetCylin
# assign material to each body in the cylinder
for body in cylinder:
body.bodyMat = material
# add cylinder to simulation
O.bodies.
# Define cylinder with funnel parameters
center1 = (0,0,height/2)
dBunker = 0.4
dOutput = 0.102
hBunker = 0
hOutput = 0.15
hPipe = 0
# create funnel as a bunker with diameter 0.102 m, height 0.064 m
funnel = geom.facetBunke
# assign material to each body in the funnel
for body in funnel:
body.bodyMat = material
# add funnel to simulation
O.bodies.
# define sphere parameters and number of spheres
rMean1 = (0.0125+0.019)/4
rRelFuzz1 = ((0.019-
num1 = 17
rMean2 = (0.0095+0.0125)/4
rRelFuzz2 = ((0.0125-
num2 = 53
rMean3 = (0.00475+0.0095)/4
rRelFuzz3 = ((0.0095-
num3 = 1267
rMean4 = (0.00236+0.00475)/4
rRelFuzz4 = ((0.00475-
num4 = 11624
#create empty sphere packing
sp = pack.SpherePack()
# generate randomly sphere
sp.makeCloud(
sp.makeCloud(
sp.makeCloud(
sp.makeCloud(
# add the sphere pack to the simulation
sp.toSimulation
for body in O.bodies:
if not isinstance(
continue
if body.shape.radius >= rMean1 :
if body.shape.radius <= rMean1 and body.shape.radius > rMean2:
if body.shape.radius <= rMean2 and body.shape.radius > rMean3:
if body.shape.radius <= rMean3 and body.shape.radius > rMean4:
if body.shape.radius <= rMean4 :
O.engines = [
# sphere, facet, wall
# the loading plate is a wall, we need to handle sphere+sphere, sphere+facet, sphere+wall
),
# the label creates an automatic variable referring to this engine
# we use it below to change its attributes from the functions called
]
O.dt = PWaveTimeStep()
for body in O.bodies:
if not isinstance(
continue
if body.shape.radius == 0.0125/2: #SP
if body.shape.radius == 0.0095/2: #SP1
if body.shape.radius == 0.00475/2: #SP2
if body.shape.radius == 0.00236/2: #SP3
# the following checkUnbalanced, unloadPlate and stopUnloading functions are all called by the 'checker'
# (the last engine) one after another; this sequence defines progression of different stages of the
# simulation, as each of the functions, when the condition is satisfied, updates 'checker' to call
# the next function when it is run from within the simulation next time
# check whether the gravity deposition has already finished
# if so, add wall on the top of the packing and start the oedometric test
def checkUnbalanced():
# at the very start, unbalanced force can be low as there is only few contacts, but it does not mean the packing is stable
if O.iter < 15000:
return
# add plate at the position on the top of the packing
# the maximum finds the z-coordinate of the top of the topmost particle
O.bodies.
global plate # without this line, the plate variable would only exist inside this function
plate = O.bodies[-1] # the last particles is the plate
# Wall objects are "fixed" by default, i.e. not subject to forces
# prescribing a velocity will therefore make it move at constant velocity (downwards)
plate.state.vel = (0, 0, -.8)
# start plotting the data now, it was not interesting before
O.engines = O.engines + [PyRunner(
# next time, do not call this function anymore, but the next one (unloadPlate) instead
checker.command = 'unloadPlate()'
def unloadPlate():
# if the force on plate exceeds maximum load, start unloading
if abs(O.forces.
plate.state.vel *= -1
# next time, do not call this function anymore, but the next one (stopUnloading) instead
checker.command = 'stopUnloading()'
def stopUnloading():
if abs(O.forces.
# calculate the volume of the packing
num_spheres = 0
for b in O.bodies:
if isinstance(b.shape, yade.wrapper.
# print the number of spheres and volume of packking
print("V Packing:", "{:e}".
# print the height of the plate
top = abs(plate.
# calculate the volume of the cylinder
# calculate the porosity and porosity percentage
print("V Cylinder:", "{:e}".
# O.tags can be used to retrieve unique identifiers of the simulation
# if running in batch, subsequent simulation would overwrite each other's output files otherwise
# d (or description) is simulation description (composed of parameter values)
# while the id is composed of time and process number
O.pause()
def addPlotData():
if not isinstance(
plot.addData()
return
Fz = O.forces.
plot.addData(
# besides unbalanced force evolution, also plot the displacement-force diagram
plot.plots = {'i': ('unbalanced',), 'w': ('Fz',)}
plot.plot()
Is there any mistake is the porosity calculation process?
Question information
- Language:
- English Edit question
- Status:
- Solved
- For:
- Yade Edit question
- Assignee:
- No assignee Edit question
- Solved by:
- Karol Brzezinski
- Solved:
- Last query:
- Last reply: