New question #701174 on Yade:
https://answers.launchpad.net/yade/+question/701174
this is an example of a triaxial shear test
from yade import pack,plot,qt,export
IsoSigma = -100.
O.periodic=True
# define material
idSand=O.materials.append(FrictMat(young=100e3,poisson=1.,frictionAngle=radians(20),density=2650,label='sand'))
# create particles
sp1=pack.SpherePack()
sp1.makeCloud(maxCorner=(0.02, 0.02, 0.02), psdSizes=[0.00017, 0.000191,
0.0002285, 0.00026, 0.000292, 0.000325, 0.00035], psdCumm=[0.0, 0.1, 0.3, 0.5,
0.6, 0.9, 1], periodic=True,num=5000,seed=1)
sp1.toSimulation(color=(0,0,1),material=idSand)
show how to use makeClumpTemplate():
#dyad:
relRadList1 = [1,1]
relPosList1 = [[0.4,0,0],[-0.4,0,0]]
#peanut:
#relRadList2 = [.5,1,.5]
#relPosList2 = [[1,0,0],[0,0,0],[-1,0,0]]
#stick:
#relRadList3 = [1,1,1,1,1]
#relPosList3 = [[0,1,0],[0,2,0],[0,3,0],[0,4,0],[0,5,0]]
templates= []
templates.append(clumpTemplate(relRadii=relRadList1,relPositions=relPosList1))
#templates.append(clumpTemplate(relRadii=relRadList2,relPositions=relPosList2))
#templates.append(clumpTemplate(relRadii=relRadList3,relPositions=relPosList3))
show how to use replaceByClumps():
#replace by 50% dyads, 30% peanuts and 10% sticks:
O.bodies.replaceByClumps(templates,[1.],discretization=10)
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]),
# PyRunner(command='fabric()',iterPeriod=1),
GlobalStiffnessTimeStepper(),
NewtonIntegrator(damping=0.2),
PeriTriaxController(goal=(IsoSigma,IsoSigma,IsoSigma), # Vector6 of
prescribed final values
stressMask=7,
dynCell=True,
maxStrainRate=(3.e+0,3.e+0,3.e+0),
maxUnbalanced=0.0001,
relStressTol=1.e-3,
doneHook='Finished()',
label='p3d'
),
# VTKRecorder(fileName='3d-vtk-',recorders=['all'],iterPeriod=1),
# PyRunner(command='strains()',iterPeriod=1),
# PyRunner(command='balls()',iterPeriod=1),
#PyRunner(command='plotAddData()',iterPeriod=100),
]
'''
def plotAddData():
plot.addData(
iter=O.iter,iter_=O.iter,
sxx=p3d.stress[0],syy=p3d.stress[1],szz=p3d.stress[2],
exx=O.cell.size[0],eyy=O.cell.size[1],ezz=O.cell.size[2],
Z=avgNumInteractions(),
Zm=avgNumInteractions(skipFree=True),
poros=porosity(),
# poros=voxelPorosity(500,(0,0,0),O.cell.size),
unbalanced=utils.unbalancedForce(),
t=O.time,
gWork=O.energy['gravWork'],
Ep=O.energy['elastPotential'],
Edamp=O.energy['nonviscDamp'],
Ediss=O.energy['plastDissip'],
Ekin=utils.kineticEnergy(),
Etot=O.energy.total(),**O.energy
)
plot.saveDataTxt('macroFile',vars=('t','exx','eyy','ezz','sxx','syy','szz','Z','Zm','poros'))
plot.saveDataTxt('energyFile',vars=('t','Etot','unbalanced','gWork','Edamp','Ekin'))
O.trackEnergy=True
# plotting
plot.live=True
plot.plots={'iter':('sxx','syy','szz'),'iter_':('exx','eyy','ezz'), '
iter':('unbalanced'),
#energy plot
' iter ':(O.energy.keys,None,'Etot')
}
def fabric():
fileNomb='fabric'+str(O.iter)
f = open(fileNomb,'w')
f.write('id1 id2 x_cp y_cp z_cp n_x n_y n_z Fn_x Fn_y Fn_z Fs_x Fs_y
Fs_z ovp\n')
for i in O.interactions:
if not i.isReal: continue
tacts = i.geom.contactPoint
nmls = i.geom.normal
ovp = i.geom.penetrationDepth
Fn = i.phys.normalForce
Fs = i.phys.shearForce
f.write('%-16d %-16d %-16g %-16g %-16g %-16g %-16g %-16g %-16g
%-16g %-16g %-16g %-16g %-16g %-16g\n'%(i.id1,
i.id2,tacts[0],tacts[1],tacts[2],nmls[0],nmls[1],nmls[2],Fn[0],Fn[1],Fn[2],Fs[0],Fs[1],Fs[2],ovp))
f.close
def balls():
fileCalled='ball'+str(O.iter)
export.text(fileCalled)
def strains():
fileNome='strain'+str(O.iter)
f = open(fileNome,'w')
f.write('b.id velx vely velz rotx roty rotz\n')
for b in O.bodies:
if isinstance(b.shape,Sphere):
vels = b.state.vel
rots = b.state.rot()
f.write('%-16d %-16g %-16g %-16g %-16g %-16g
%-16g\n'%(b.id, vels[0],vels[1],vels[2],
rots[0],rots[1],rots[2]))
f.close
'''
def Finished():
