yade-users team mailing list archive

[JIPEIQI <question668191@xxxxxxxxxxxxxxxxxxxxx>]

Question #668191 on Yade changed:
https://answers.launchpad.net/yade/+question/668191

    Status: Solved => Open

JIPEIQI is still having a problem:
Hi Jan, I found the makeCloud function seems doesn't create non-
overlapping cloud. The script I use is followed below:


#############################
from yade import ymport,pack,export
import sys
#surf=ymport.stl('disc.stl',fixed=False)
#O.bodies.append(surf)
key='a'
mv_x=200e-3
mv_y=200e-3
mv_z=200e-3
sxx=-2e5
targetPorosity = 0.50
rlow=3e-3
rhi=rlow*1.5
rmean=(rlow+rhi)/2.0
rfuzz=(rhi-rlow)/2.0/rmean
rmean2=rmean/2.0
numBall=int(mv_x*mv_y*mv_z/(4.0/3.0*pi*rmean**3)*(1-targetPorosity))
print numBall
O.materials.append(FrictMat(young=20e9,poisson=0.5,frictionAngle=0,density=2600,label='spheres'))
O.materials.append(FrictMat(young=20e9,poisson=0.5,frictionAngle=0,density=2600,label='walls'))


mn,mx=Vector3(-mv_x/2.0,-mv_y/2.0,mv_z),Vector3(mv_x/2.0,mv_y/2.0,0)
walls=aabbWalls([mn,mx],thickness=0,material="walls")
wallIds=O.bodies.append(walls)
center=Vector3(0,0,mv_z/2.0)

sp=pack.SpherePack()
sp.makeCloud(mn,mx,rmean2,rfuzz,numBall,False, 0.95,seed=1)
sp.toSimulation(material="spheres")
for b in O.bodies:
	if isinstance(b.shape,Sphere):
		b.state.pos=(b.state.pos-center)*1+center
export.text('particles.txt',-1)

def getPoro():
	Vol=0.0
	for o in O.bodies:
		if isinstance(o.shape,Sphere):
			Vol+=o.shape.radius**3*4.0/3.0*pi
	return 1-Vol/mv_x/mv_y/mv_z
radiusMult = ( (1.0 -targetPorosity ) / (1.0 - getPoro()) )**(1.0/3.0)

print radiusMult

def expand(mult):
	print "expanding radius * %f" % mult
	for b in O.bodies:
		if isinstance(b.shape,Sphere):
			b.shape.radius*=mult
			b.state.mass=2600*pi*4.0/3.0*b.shape.radius**3
def calm():
	for b in O.bodies:
		if isinstance(b.shape,Sphere):
			b.state.vel=Vector3(0,0,0)
			b.state.angVel=Vector3(0,0,0)
print getPoro()
triax=ThreeDTriaxialEngine(
	maxMultiplier=1.+2.0e4/20e9,
	finalMaxMultiplier=1.+2.0e3/20e9,
	thickness = 0,
	stressControl_1 = True,
	stressControl_2 = True,
	stressControl_3 = True,
	## Independant stress values for anisotropic loadings
	goal1=sxx,
	goal2=sxx,
	goal3=sxx,
	internalCompaction=True,
	Key=key,
)
newton=NewtonIntegrator(damping=0.7)

O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Wall_Aabb()]),
	InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
		[Ip2_FrictMat_FrictMat_FrictPhys()],
		[Law2_ScGeom_FrictPhys_CundallStrack()]
	),
	newton,
	VTKRecorder(fileName='vtk/gravity-',recorders=['all'],label='vvtk',iterPeriod=1000)
]
O.dt=.5*PWaveTimeStep()
def equ(n):
	M=radiusMult**(1/float(n))
	for nn in range(0,n):
		expand(M)
		print 'step1: current poro: ', getPoro()
		for nn in range(0,50):
			O.run(5,True)
			calm()
		while 1:
			O.run(1000, True)
  	#the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium
			unb=unbalancedForce()
			print 'unbalanced force:', unb
			calm
			if unb<1e-4:
				break
##################################################


#for ip in O.interactions:
#	print O.bodies[ip.id1].state.pos
#	print O.bodies[ip.id1].state.mass
#	print O.bodies[ip.id2].state.pos
#	print O.bodies[ip.id2].state.mass
#	print ip.phys.kn
#	print ip.phys.normalForce
#	break
yade.qt.Controller(), yade.qt.View()
#equ(5)
O.run(1,True)
def hh():
	for ip in O.interactions:
		#print O.bodies[ip.id1].state.vel
		#print O.bodies[ip.id2].state.vel
		print ip.geom.penetrationDepth/rmean2
hh()
################################################################

the print output is usually 0.5-0.9, which is very large to make the
particles explodes out of the box

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.