New question #706667 on Yade:
https://answers.launchpad.net/yade/+question/706667
I use the code by Luc. But there is an error. How to solve?
from yade import ymport, utils , plot
#---------------- DEFINITION OF SIMULATION'S PARAMETERS
#### packing (previously constructed)
PACKING='121_1k.spheres'
#### Boundary Conditions
confinement=-1e6
strainRate=-0.02
#### name of output files
OUT=PACKING+'_1MPa_r0.02'
#### Simulation Control
saveData=10 # data record interval
iterMax=50000 # maximum number of iterations
saveVTK=iterMax/5 # Vtk files record interval
#### Material microproperties -> Lac du Bonnet granite (cf. A DEM model for
soft and hard rock, Scholtes & Donze, JMPS 2013)
intR=1.5 # allows near neighbour interaction and defines coordination number
K=13 (needs to be adjusted for every packing -> preprocessing needed)
DENS=4000 # could be adapted to match material density:
dens_DEM=dens_rock*(V_rock/V_particles)=dens_rock*1/(1-poro_DEM) -> poro?
YOUNG=68e9
FRICT=10
ALPHA=0.4
TENS=8e6
COH=160e6
#### material definition
def sphereMat(): return JCFpmMat(type=1,density=DENS,young=YOUNG,poisson =
ALPHA,frictionAngle=radians(FRICT),tensileStrength=TENS,cohesion=COH)
def wallMat(): return
JCFpmMat(type=0,density=DENS,young=YOUNG,frictionAngle=radians(0))
#### preprocessing to get dimensions
O.bodies.append(ymport.text(PACKING,scale=1.,shift=Vector3(0,0,0),material=sphereMat))
dim=utils.aabbExtrema()
xinf=dim[0][0]
xsup=dim[1][0]
X=xsup-xinf
yinf=dim[0][1]
ysup=dim[1][1]
Y=ysup-yinf
zinf=dim[0][2]
zsup=dim[1][2]
Z=zsup-zinf
R=0
Rmax=0
numSpheres=0.
for o in O.bodies:
if isinstance(o.shape,Sphere):
numSpheres+=1
R+=o.shape.radius
Rmean=R/numSpheres
#### IMPORTANT LINE HERE
O.reset() # all previous lines were for getting dimensions of the packing to
create walls at the right positions (below) because walls have to be genrated
after spheres for FlowEngine
#### now we construct the scene with right dimensions (because walls have to be
imported before spheres for certain engines)
### walls
mn,mx=Vector3(xinf+0.1*Rmean,yinf+0.1*Rmean,zinf+0.1*Rmean),Vector3(xsup-0.1*Rmean,ysup-0.1*Rmean,zsup-0.1*Rmean)
walls=utils.aabbWalls(oversizeFactor=1.5,extrema=(mn,mx),thickness=min(X,Y,Z)/100.,material=wallMat)
wallIds=O.bodies.append(walls)
### packing
O.bodies.append(ymport.text(PACKING,scale=1.,shift=Vector3(0,0,0),material=sphereMat))
#---------------- ENGINES ARE DEFINED HERE
#### triaxial Engine
triax=TriaxialStressController(
internalCompaction=False
)
#### simulation is defined here (DEM loop, interaction law, servo control,
recording, etc...)
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb(aabbEnlargeFactor=intR,label='Saabb')]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intR,label='SSgeom'),Ig2_Box_Sphere_ScGeom()],
[Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
[Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,Key=OUT,label='interactionLaw')]
),
triax,
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.5,defaultDt=0.1*utils.PWaveTimeStep()),
NewtonIntegrator(damping=0.5,label="newton"),
PyRunner(iterPeriod=int(saveData),initRun=True,command='recorder()',label='data'),
VTKRecorder(iterPeriod=int(saveVTK),initRun=True,fileName=OUT+'-',recorders=['spheres','jcfpm','cracks'],Key=OUT,label='vtk')
]
#### custom recording functions
tensCks=shearCks=0
e10=e20=e30=0
def recorder():
global tensCks,shearCks,e10,e20,e30
tensCks=0
shearCks=0
for o in O.bodies:
if isinstance(o.shape,Sphere):
tensCks+=o.state.tensBreak
shearCks+=o.state.shearBreak
yade.plot.addData(
t=O.time,i=O.iter,e1=triax.strain[0]-e10,e2=triax.strain[1]-e20,e3=triax.strain[2]-e30,s1=triax.stress(triax.wall_right_id)[0],s2=triax.stress(triax.wall_top_id)[1],s3=triax.stress(triax.wall_front_id)[2],tc=0.5*tensCks,sc=0.5*shearCks,unbF=utils.unbalancedForce())
plot.saveDataTxt(OUT)
# if you want to plot during simulation
plot.plots={'i':('s1','s2','s3')}
plot.plot()
#---------------- SIMULATION STARTS HERE
#### manage interaction detection factor during the first timestep and then set
default interaction range (intRadius=1)
#O.step()
### initializes the interaction detection factor
SSgeom.interactionDetectionFactor=-1.
Saabb.aabbEnlargeFactor=-1.
#### APPLYING ISOTROPIC LOADING
triax.stressMask=7
triax.goal1=confinement
triax.goal2=confinement
triax.goal3=confinement
triax.max_vel=0.001
#while 1:
#if confinement==0:
#O.run(1000,True) # to stabilize the system
#break
#O.run(100,True)
#unb=unbalancedForce()
#meanS=abs(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3
#print ('unbalanced force:',unb,' mean stress: ',meanS)
#if unb<0.005 and abs(meanS-abs(confinement))/abs(confinement)<0.001:
#O.run(1000,True) # to stabilize the system
#e10=triax.strain[0]
#e20=triax.strain[1]
#e30=triax.strain[2]
#break
#### APPLYING DEVIATORIC LOADING ALONG Y AXIS
#triax.stressMask=5
#triax.goal1=confinement
#triax.goal2=strainRate
#triax.goal3=confinement
#triax.max_vel=1
#O.run(iterMax)
--
You received this question notification because your team yade-users is
an answer contact for Yade.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
