yade-users team mailing list archive

[JIPEIQI <question668191@xxxxxxxxxxxxxxxxxxxxx>]

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

Hello everyone, I 'm using script like [1] to generate a dense sample which will be later cemented as intact rock.  The script is like 

from yade import ymport,pack,export
import sys
#surf=ymport.stl('disc.stl',fixed=False)
#O.bodies.append(surf)
key='test'
mv_x=200e-3
mv_y=200e-3
mv_z=200e-3
targetPorosity = 0.40
rlow=3e-3
rhi=rlow*1.5
rmean=(rlow+rhi)/2.0
rfuzz=(rhi-rlow)/2.0/rmean
rmean2=rmean/1.1
numBall=int(mv_x*mv_y*mv_z/(4.0/3.0*pi*rmean**3)*(1-targetPorosity))
print numBall
#O.materials.append(FrictMat(young=1e9,poisson=0.5,frictionAngle=0,density=2600,label='spheres'))
O.materials.append(FrictMat(young=1e9,poisson=0.5,frictionAngle=0,density=0,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)
try:
	O.bodies.append(ymport.text('particles.txt'))
	print len(O.bodies),numBall
	a=raw_input('continue y/n')
	if a=='n':
		sys.exit(0)
except:
	sp=pack.SpherePack()
	sp.makeCloud(mn,mx,rmean2,rfuzz,numBall,False, 0.95,seed=1)
	sp.toSimulation()
	export.text('particles.txt',-1)
	print 'fooled'

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():
	print "expanding radius * %f" % radiusMult
	for b in O.bodies:
		if isinstance(b.shape,Sphere):
			b.shape.radius*=radiusMult
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.005,
	finalMaxMultiplier=1.002,
	thickness = 0,
	stressControl_1 = True,
	stressControl_2 = True,
	stressControl_3 = True,
	## Independant stress values for anisotropic loadings
	goal1=-1e6,
	goal2=-1e6,
	goal3=-1e6,
	internalCompaction=True,
	Key=key,
)
newton=NewtonIntegrator(damping=0.7)

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()]
	),
	triax,
	newton,
	VTKRecorder(fileName='vtk/gravity-',recorders=['all'],label='vvtk',iterPeriod=400)
]
O.dt=.5*PWaveTimeStep()
print 'current poro: ', getPoro()
#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,' mean stress: ',triax.meanStress
	#print 'current poro: ', getPoro()
	#if unb<1e-4 and abs(-10000-triax.meanStress)/10000<0.001:
		#break
O.run(1.True)
O.save('disc.yade')
######################################################################

The particles will explodes with nothing left as I can see nothing in the paraview or in the plot window
[1]https://github.com/yade/trunk/tree/master/examples/triax-tutorial

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