yade-users team mailing list archive

[Karol Brzezinski <question706766@xxxxxxxxxxxxxxxxxxxxx>]

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

Karol Brzezinski proposed the following answer:
Hi Leonard,

Try this (please note that I am using python 3, so the call pf print function is modified):
##########
################### 2D MWE #############
from __future__ import division
from yade import pack, plot, export, ymport

num_spheres=1000

rate=-0.01
damp=0.6
stabilityThreshold=0.001
young=5e6
confinement=6.7e3

mn,mx=Vector3(0,0,0.02),Vector3(1,2,0.02)

O.materials.append(FrictMat(young=young, poisson=0.5, frictionAngle=radians(30), density=2600, label='spheres'))
O.materials.append(FrictMat(young=young, poisson=0.5, frictionAngle=0, density=0, label='walls'))
walls=aabbWalls([Vector3(0,0,0),Vector3(1,2,0.04)],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

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

O.bodies.append([sphere(center,rad,material='spheres') for center,rad in
sp])

Gl1_Sphere.quality=3

for b in O.bodies:
        if isinstance(b.shape,Sphere):
                b.state.blockedDOFs='zXY'
                b.shape.color=[1,1,1]

triax=TriaxialStressController(
        thickness = 0,
        stressMask = 7,
        internalCompaction=False,
)

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,
        newton
]

triax.goal1=triax.goal2=triax.goal3=-confinement

while 1:
  O.run(1000, True)
  unb=unbalancedForce()
  print('unbF:',unb,' meanStress: ',-triax.meanStress,'top:',-triax.stress(triax.wall_top_id)[1],'left:',-triax.stress(triax.wall_left_id)[0],'front:',-triax.stress(triax.wall_front_id)[2])
  if unb<stabilityThreshold and abs(-confinement-triax.meanStress)/confinement<0.0001:
          break

print("### state 1 completed ###")


######## export files and run simulation 
export.text('state_0.txt')## export files state 0
O.run(500,True)
export.text('state_1.txt')## export files state 1


###################### 
layerThickness = 0.015*2 # this is more or less diameter of your particle
O.bodies.clear()
for i in range(3):
    bb = ymport.text('state_0.txt', shift = Vector3(0,0,i*layerThickness))
    O.bodies.append(bb)
TW=TesselationWrapper()
TW.triangulate()
TW.computeVolumes()
TW.setState(0)

O.bodies.clear()
for i in range(3):
    bb = ymport.text('state_1.txt', shift = Vector3(0,0,i*layerThickness))
    O.bodies.append(bb)
TW.setState(1)
TW.defToVtk("strain_pseudo2D.vtk")

#######
Cheers,
Karol

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