yade-users team mailing list archive

[Lukas Wasner <question703081@xxxxxxxxxxxxxxxxxxxxx>]

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.