yade-users team mailing list archive

[kelaogui <question216150@xxxxxxxxxxxxxxxxxxxxx>]

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

kelaogui posted a new comment:
Yeah,
    Thank you for your suggestion. Here is my code. 


from yade import pack
from yade import *

num_spheres=2000
## corners of the initial packing
mn,mx=Vector3(0,0,0),Vector3(1,1,1)
thick = 0.01
compFricDegree = 2
rate=0.2
damp=0.1
stabilityThreshold=0.001
key='_define_a_name_'

## create material #0, which will be used as default
O.materials.append(FrictMat(young=5e6,poisson=0.5,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=5e6,poisson=0.5,frictionAngle=0,density=0,label='walls'))

## create walls around the packing
walls=utils.aabbWalls([mn,mx],thickness=thick,material='walls')
wallIds=O.bodies.append(walls)

sp=pack.SpherePack()
sp.makeCloud(mn,mx,-1,0.3333,num_spheres,False, 0.95)


volume = (mx[0]-mn[0])*(mx[1]-mn[1])*(mx[2]-mn[2])
mean_rad = pow(0.09*volume/num_spheres,0.3333)

clumps=True
if clumps:
	c1=pack.SpherePack([((-0.2*mean_rad,0,0),0.5*mean_rad),((0.2*mean_rad,0,0),0.5*mean_rad)])
	sp.makeClumpCloud((0,0,0),(1,1,1),[c1],periodic=False)
	O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in sp])
	standalone,clumps=sp.getClumps()
	for clump in clumps:
		O.bodies.clump(clump)
		for i in clump[1:]: O.bodies[i].shape.color=O.bodies[clump[0]].shape.color
	#sp.toSimulation()
else:
	O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in sp])

O.dt=.5*utils.PWaveTimeStep() # initial timestep, to not explode right away
O.usesTimeStepper=True
triax=ThreeDTriaxialEngine(
	#maxMultiplier=1.005,
	#finalMaxMultiplier=1.002,
	thickness = thick,
	stressControl_1 = True,
	stressControl_2 = True,
	stressControl_3 = True,
        strainRate1=0.01,
        strainRate2=0.01,
        strainRate3=0.01,
	## The stress used for (isotropic) internal compaction
	sigma_iso = 10000,
	## Independant stress values for anisotropic loadings
	sigma1=10000,
	sigma2=10000,
	sigma3=10000,
	internalCompaction=False,
	Key=key,
)

newton=NewtonIntegrator(damping=damp)

O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()],verletDist=-mean_rad*0.06),
	InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
		[Ip2_FrictMat_FrictMat_FrictPhys()],
		[Law2_ScGeom_FrictPhys_CundallStrack(label='cundall')]
	),
	GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8, defaultDt=4*utils.PWaveTimeStep()),
	triax,
	TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+key),
	newton
]

#Display spheres with 2 colors for seeing rotations better
Gl1_Sphere.stripes=0
yade.qt.Controller(), yade.qt.View()

while 1:
  O.run(100, True)
  #the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium
  unb=unbalancedForce()
  meanS=(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3
  sidelen1=O.bodies[1].state.se3[0][0]-O.bodies[0].state.se3[0][0]-thick
  print 'unbalanced force:',unb,' mean stress: ',meanS,'side:  ',sidelen1
  if unb<stabilityThreshold and abs(meanS-triax.sigma_iso)/triax.sigma_iso<0.001:
    print "###      Isotropic state saved      ###"
    break

-- 
You received this question notification because you are a member of
yade-users, which is an answer contact for Yade.