Question #703081 on Yade changed:
https://answers.launchpad.net/yade/+question/703081
Status: Answered => Solved
Lukas Wasner confirmed that the question is solved:
Hello,
I tried a new approach that allowed me to achieve the target porosity, at least
for a cuboid predicate.
For a cylindrical predicate, with the same approach, there is a higher
porosity at the beginning, but it slowly decreases with increasing
iteration (I did not run the simulation completely due to time
constraints).
The approach is based on Bettina's hints that spheres from the sphere
packing are inserted into the predicate; however, spheres that intersect
the predicate are completely removed. This leads to an increased
porosity, which is basically a boundary value problem.
The idea is this:
the predicate is chosen larger than the facet volume.
The predicate boundary is increased by the maximum existing sphere radius in
the sphere packing in each spatial dimension by a factor of 2.
If the predicate is larger than the actual volume, 2 things happen:
a) spheres that lie in the predicate but intersect the facet volume are not
deleted. This increases the number of spheres in the volume and thus decreases
the porosity.
b) spheres with a sphere radius smaller than the maximum sphere radius can lie
outside the facet volume (since the predicate is larger than it).
Now the following is done:
Since the spatial facet boundaries are known, all spheres are deleted from the
filtered sphere packing that have their center (state.pos) outside these
boundaries (condition b))
Spheres whose centers are within the boundaries but intersect the facets will
slip into the volume when the simulation is started (condition a)).
The start of the simulation is somewhat chaotic, but after an
equilibrium has been established, the actual simulation (e.g. oedometer
test) can be started.
--------------------------------------------------------------------------------------------------------
generating sphere pack with target porosity:
-----------------------------------------------------------------------------------------------------
readParamsFromTable(factor=5, sig3=-10000)
from yade import pack, plot
from yade.params.table import *
targetPorosity = 0.4 #the porosity we want for the packing
damp=0.2 # damping coefficient
stabilityThreshold=0.01 # we test unbalancedForce against this value in
different loops (see below)
young=5e6 # contact stiffness
mn,mx=Vector3(0,0,0),Vector3(.1,.1,.1) # corners of the initial packing
compFricDegree=20
psdSizes,psdCumm=[factor*.0003,factor*.0004,factor*.0006,factor*.0008],[0,.5,.5,1]
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)
sp=pack.SpherePack()
sp.makeCloud(mn,mx,psdSizes=psdSizes,psdCumm=psdCumm,distributeMass=True,seed=1)
sp.toSimulation(material='spheres')
print('Anzahl Kugeln:', len(sp))
triax=TriaxialStressController(
thickness = 0,
stressMask = 7,
internalCompaction=False, # If true the confining pressure is generated
by growing particles
goal1=sig3,
goal2=sig3,
goal3=sig3
)
newton=NewtonIntegrator(damping=damp)
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
),
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
triax,
#TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key),
newton,
PyRunner(command='checkUnbalanced()', iterPeriod=1000, label='checker'),
PyRunner(command='addPlotData()', iterPeriod=200),
]
def checkUnbalanced():
unb=unbalancedForce()
print('unbalanced Force: ', unb, ' mean stress: ', triax.meanStress)
if O.iter> 1000 and unb<stabilityThreshold and
abs((triax.goal1-triax.meanStress)/triax.goal1)<.001:
checker.iterPeriod=500
checker.command='reduceFric()'
def reduceFric():
if triax.porosity>targetPorosity:
compFricDegree=.95*compFricDegree
setContactFriction(radians(comFricDegree))
print('Friction: ', compFricDegree, 'Porosity: ',triax.porosity)
#sys.stdout.flush()
else:
#sp.fromSimulation()
#sp.save(f'STRESSED_Pack_factor_{factor}_porosity_{triax.porosity}.txt')
print('Destress')
checker.iterPeriod=150
checker.command='destress()'
def destress():
triax.goal1=sig3*.01
triax.goal2=sig3*.01
triax.goal3=sig3*.01
if abs((triax.goal1-triax.meanStress)/triax.goal1)<.001 and
triax.porosity<targetPorosity:
print('Porosity Reached: ',triax.porosity)
sp.fromSimulation()
#sp.save(f'DESTRESSED_Pack_factor_{factor}_porosity_{triax.porosity}.txt')
sp.save('testPack.txt')
O.pause()
def addPlotData():
plot.addData(unbalanced=unbalancedForce(), Porosity=triax.porosity,
Sigma=triax.goal1, MeanStress=triax.meanStress, i=O.iter)
plot.plots = {'i': ('unbalanced'), 'i ': ('Porosity'), 'i ': ('Sigma',
'MeanStress')}
plot.plot()
O.run()
waitIfBatch()
-----------------------------------------------------------------
filter sphere pack to prediacte
-----------------------------------------------------------------
from yade import pack, plot
import numpy as np
#from yade.params.table import *
height=.02
#radiusoed=.0356825
radiusoed=.03
overlap=True
Cyl=False
sp=SpherePack()
sp.load('testPack.txt')
maxR=max(s[1] for s in sp)
minR=min(s[1] for s in sp)
meanR=np.average([s[1] for s in sp])
if Cyl==False:
O.bodies.append(geom.facetBox(center=(radiusoed,radiusoed,height/2),extents=(radiusoed,radiusoed,height/2,),wallMask=63))
if overlap==False:
pred=pack.inAlignedBox(minAABB=(0,0,0),maxAABB=(radiusoed*2,radiusoed*2,height))
else:
pred=pack.inAlignedBox(minAABB=(0-maxR,0-maxR,0-maxR),maxAABB=(radiusoed*2+maxR,radiusoed*2+maxR,height+maxR))
else:
O.bodies.append(geom.facetCylinder(center=(0, 0, height/2),
radius=radiusoed, height=height,segmentsNumber=2000, wallMask=7)) #oder6
if overlap==False:
pred=pack.inCylinder(centerBottom=(0,0,0),
centerTop=(0,0,height), radius=radiusoed)
else:
pred=pack.inCylinder(centerBottom=(0,0,0-maxR/4),
centerTop=(0,0,height+maxR/4), radius=radiusoed+maxR/3)
sp=pack.filterSpherePack(predicate=pred,spherePack=sp,returnSpherePack=True)
sp.toSimulation()
for b in O.bodies:
if isinstance(b.shape,Sphere):
b.shape.color=scalarOnColorScale(x=b.shape.radius, xmin=minR,
xmax=maxR)
print('Anzahl Körper vor löschen: ',len(O.bodies))
ListKugeln=[]
for b in O.bodies:
if isinstance(b.shape,Sphere):
ListKugeln.append(b.id)
ListDel=[]
if Cyl==False:
for b in O.bodies:
if isinstance(b.shape, Sphere):
if b.state.pos[0]<0:
ListDel.append(b.id)
if b.state.pos[0]>radiusoed*2:
ListDel.append(b.id)
if b.state.pos[1]<0:
ListDel.append(b.id)
if b.state.pos[1]>radiusoed*2:
ListDel.append(b.id)
if b.state.pos[2]<0:
ListDel.append(b.id)
if b.state.pos[2]>height:
ListDel.append(b.id)
else:
for b in O.bodies:
if isinstance(b.shape, Sphere):
if b.state.pos[2]<0:
ListDel.append(b.id)
if b.state.pos[2]>height:
ListDel.append(b.id)
if sqrt(b.state.pos[0]**2+b.state.pos[1]**2)>radiusoed:
ListDel.append(b.id)
ListKugeln=list(set(ListKugeln))
ListDel=list(set(ListDel))
for x in ListDel:
O.bodies.erase(x)
print('Anzahl aller Kugeln: ', len(ListKugeln))
print('Anzahl gelöschter Kugeln: ',len(ListDel))
print('Anzahl alle Kugeln nach Bereinigung: ', len(ListKugeln)-len(ListDel))
print('Porosity after filter: ', porosity())
O.engines = [
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb(),
Bo1_Wall_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom(),
Ig2_Box_Sphere_ScGeom(), Ig2_Wall_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
),
NewtonIntegrator(gravity=(0, 0, -9.81), damping=0.9),
# the label creates an automatic variable referring to this engine
# we use it below to change its attributes from the functions called
PyRunner(command='checkUnbalanced()', iterPeriod=20, label='checker'),
PyRunner(command='addPlotData()', iterPeriod=20),
]
O.dt =PWaveTimeStep()
def checkUnbalanced():
if unbalancedForce()<.005:
print('Porenanteil n : ',porosity())
print('Porenzahl e: ', porosity()/(1-porosity()))
ListLast=[]
for b in O.bodies:
if isinstance(b.shape, Sphere):
ListLast.append(b.id)
print('Anzahl Kugeln nach Bewegung: ',len(ListLast))
print('Kugeln die verloren gegangen sind:
',len(ListKugeln)-len(ListDel)-len(ListLast))
O.pause()
def addPlotData():
plot.addData(unbalanced=unbalancedForce(),Porosity=porosity(),i=O.iter)
plot.plots = {'i': ('unbalanced'), 'i ': ('Porosity')}
plot.plot()
-------------------------------------------------------------------
Note that at the beginning of the filter predicate script there are boolean
variables that can be used to set the predicate form (Cyl=False for cuboid,
Cyl=True for cylinder).
In addition, the particle-facet overlap can also be switched on and off here
(overlap=True).
best regards,
Lukas
--
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 : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
