yade-users team mailing list archive

[Luc Scholtès <question556907@xxxxxxxxxxxxxxxxxxxxx>]

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

Luc Scholtès proposed the following answer:
Hi Weijy,

Sorry for not working directly on your script but I found it easier to
actually send you a script that I used to obtain the results presented
in the paper you mentioned in your first post.

As I told you, I usually work with predefined dense packings (I use a
dedicated script for that). First because I generate them using the same
procedure and make thus sure they have similar statistical properties
(important for the near neighbour intR feature). Second because when I
run several triaxial tests, I prefer to use always the same packing
whatever the confining stress to avoid any bias related to the packing.

Please find below the script (you'll have to adjust the intR parameter
to the packing, I just defined a default value here).

I hope it will help you to get what you want.

Luc

---

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:
	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.