yade-users team mailing list archive

[Luc Scholtès <question661250@xxxxxxxxxxxxxxxxxxxxx>]

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

    Status: Open => Answered

Luc Scholtès proposed the following answer:
Wow... so many lines.

I would suggest that you use the following first batch of lines to get
used to yade first. If you copy paste these lines in a file and name it
"uniaxialCompression.py", you should then be able to run it with "yade
uniaxialCompression.py". Basically, it simulates a uniaxial compression
with the JCFPM model. Then, if you want to plot some data from the
simulation, you can copy the second batch of lines, copy paste them in a
file name "plotUniaxialCompression.py" and then launch it with "python
plotUniaxialCompression.py".

Then, once you know how to run and exploit this simulation, you'll be
able to modify the script with the contact law that suits you the most.
If not, please ask a simple question with minimum lines and clear
(short) title.

Luc

--- uniaxial compression:

# -*- coding: utf-8 -*-
from yade import pack, plot

################# SIMULATIONS DEFINED HERE

#### packing (previously constructed)
OUT='compressionTest_JCFPM'

#### Simulation Control
rate=-0.01 #deformation rate 
iterMax=10000 # maximum number of iterations 
saveVTK=2000 # saving output files for paraview

#### Material microproperties 
intR=1.1 # allows near neighbour interaction (can be adjusted for every packing)
DENS=2500 # could be adapted to match material density: dens_DEM=dens_rock*(V_rock/V_particles)=dens_rock*1/(1-poro_DEM) -> packing porosity as to be computed? 
YOUNG=20e9 
FRICT=7
ALPHA=0.1
TENS=1e6 
COH=1e6

#### material definition
def sphereMat(): return JCFpmMat(type=1,density=DENS,young=YOUNG,poisson=ALPHA,frictionAngle=radians(FRICT),tensileStrength=TENS,cohesion=COH)

#### create the specimen
pred=pack.inCylinder((0,0,0),(0,1,0),0.25)
O.bodies.append(pack.regularHexa(pred,radius=0.025,gap=0.,material=sphereMat))

R=0
Rmax=0
nbSpheres=0.
for o in O.bodies:
 if isinstance(o.shape,Sphere):
   nbSpheres+=1
   R+=o.shape.radius
   if o.shape.radius>Rmax:
     Rmax=o.shape.radius
Rmean=R/nbSpheres

print 'nbSpheres=',nbSpheres,' | Rmean=',Rmean

#### boundary condition (see utils.uniaxialTestFeatures
bb=utils.uniaxialTestFeatures()
negIds,posIds,longerAxis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']
 
################# ENGINES DEFINED HERE

O.engines=[
	ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intR,label='Saabb')]),
	InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intR,label='SSgeom')],
		[Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
		[Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,Key=OUT,label='interactionLaw')]
	),
	UniaxialStrainer(strainRate=rate,axis=longerAxis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=1,blockRotations=1,setSpeeds=0,stopStrain=0.1,dead=1,label='strainer'),
	GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.5, defaultDt=utils.PWaveTimeStep()),
	NewtonIntegrator(damping=0.4,label='newton'),
	PyRunner(iterPeriod=int(100),initRun=True,command='recorder()',label='data'),
        VTKRecorder(iterPeriod=int(saveVTK),initRun=True,fileName=OUT+'-',recorders=['spheres','jcfpm','cracks'],Key=OUT,label='vtk')
]

################# RECORDER DEFINED HERE

def recorder():
    yade.plot.addData({'i':O.iter,
		       'eps':strainer.strain,
		       'sigma':strainer.avgStress,
		       'tc':interactionLaw.nbTensCracks,
		       'sc':interactionLaw.nbShearCracks,
		       'te':interactionLaw.totalTensCracksE,
		       'se':interactionLaw.totalShearCracksE,
		       'unbF':utils.unbalancedForce()})
    plot.saveDataTxt(OUT)
    
# if you want to plot during simulation
plot.plots={'i':('sigma')}
#plot.plot()

################# PREPROCESSING

#### manage interaction detection factor during the first timestep and then set default interaction range ((cf. A DEM model for soft and hard rock, Scholtes & Donze, JMPS 2013))
O.step();
### initializes the interaction detection factor
SSgeom.interactionDetectionFactor=-1.
Saabb.aabbEnlargeFactor=-1.

#### coordination number verification
numSSlinks=0
numCohesivelinks=0
for i in O.interactions:
    if isinstance(O.bodies[i.id1].shape,Sphere) and isinstance(O.bodies[i.id2].shape,Sphere):
      numSSlinks+=1
    if i.phys.isCohesive :
      numCohesivelinks+=1

print "nblinks=", numSSlinks, " | nbCohesivelinks=", numCohesivelinks, " || Kcohesive=", 2.0*numCohesivelinks/nbSpheres
  
################# SIMULATION REALLY STARTS HERE
strainer.dead=0
O.run(iterMax)

--- plot uniaxial test data:

# -*- coding: utf-8 -*-
from pylab import *

### processing function
def store(var,textFile):
    data=loadtxt(textFile,skiprows=1)
    it=[]
    eps=[]
    sig=[]
    tc=[]
    sc=[]
    te=[]
    se=[]
    ubf=[]
    for i in range(0,len(data)):
      it.append(float(data[i,1]))
      eps.append(abs(float(data[i,0])))
      sig.append(abs(float(data[i,4])))
      tc.append(float(data[i,5]))
      sc.append(float(data[i,2]))
      te.append(float(data[i,6]))
      se.append(float(data[i,3]))
      ubf.append(float(data[i,7]))
    var.append(it)
    var.append(eps)
    var.append(sig)
    var.append(tc)
    var.append(sc)
    var.append(te)
    var.append(se)
    var.append(ubf)
   
### data input
dataFile1='compressionTest_JCFPM'

a1=[]
store(a1,dataFile1)

rcParams.update({'legend.numpoints':1,'font.size': 20,'axes.labelsize':25,'xtick.major.pad':10,'ytick.major.pad':10,'legend.fontsize':20})
lw=2
ms=10

### plots

figure(0,figsize=(10,10))
ax1=subplot(1,1,1)
xlabel('iteration [-]')
ax1.plot(a1[0],[x/1e6 for x in a1[2]],'-k',linewidth=lw)
ylabel(r'$\sigma_1$ [MPa]')
ax2 = ax1.twinx()
ax2.plot(a1[0],a1[5],'-r',linewidth=lw)
ylabel('unbForce [-]')

figure(1,figsize=(10,10))
ax1=subplot(1,1,1)
xlabel(r'$\epsilon_1$ [millistrain]')
plot([x*1e3 for x in a1[1]],[x/1e6 for x in a1[2]],'-k',linewidth=lw)
ylabel(r'$\sigma_1$ [MPa]')

figure(2,figsize=(10,10))
ax1=subplot(1,1,1)
xlabel(r'$\epsilon_1$ [millistrain]')
ax1.plot([x*1e3 for x in a1[1]],[x/1e6 for x in a1[2]],'-k',linewidth=lw)
ylabel(r'$\sigma_1$ [MPa]')
ax2 = ax1.twinx()
ax2.plot([x*1e3 for x in a1[1]],a1[3],'-b',linewidth=lw)
ax2.plot([x*1e3 for x in a1[1]],a1[4],'-r',linewidth=lw)
ylabel('cumulative number of microcracks [-]')
legend(('tensile','shear'),loc=2)

figure(3,figsize=(10,10))
ax1=subplot(1,1,1)
xlabel(r'$\epsilon_1$ [millistrain]')
ax1.plot([x*1e3 for x in a1[1]],[x/1e6 for x in a1[2]],'-k',linewidth=lw)
ylabel(r'$\sigma_1$ [MPa]')
ax2 = ax1.twinx()
ax2.plot([x*1e3 for x in a1[1]],[x+y for x,y in zip(a1[5],a1[6])],'-r',linewidth=lw)
ylabel('released energy [J]')

### show or save
show()

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