Question #703169 on Yade changed:
https://answers.launchpad.net/yade/+question/703169
Status: Needs information => Open
Xu gave more information on the question:
Hi Karol,
Thank you for your reply.
This is the script i use:
######Consolidation#######
from yade import pack,plot,qt,export,ymport,utils
import numpy as np
IsoSigma = -1.e5
O.periodic=True
O.trackEnergy=True
idSand=O.materials.append(FrictMat(young=2e8,poisson=.8,frictionAngle=0.09,density=2650000,label='spheres'))
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)
v=utils.getSpheresVolume()
print utils.getSpheresVolume()
fout = file('Packing_volume.txt','w')
fout.write(str(v))
fout.close()
#### show how to use makeClumpTemplate():
relRadList2 = [1,1]
relPosList2 = [[0.4,0,0],[-0.4,0,0]]
templates= []
templates.append(clumpTemplate(relRadii=relRadList2,relPositions=relPosList2))
O.bodies.replaceByClumps(templates,[1.],discretization=10)
law=Law2_ScGeom_FrictPhys_CundallStrack(traceEnergy=True)
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[law]),
# PyRunner(command='fabric()',iterPeriod=10000),
GlobalStiffnessTimeStepper(),
NewtonIntegrator(damping=0.20),
PeriTriaxController( goal=(IsoSigma,IsoSigma,IsoSigma), # Vector6 of
prescribed final values
stressMask=7,
dynCell=True,
maxStrainRate=(5.e+0,5.e+0,5.e+0),
maxUnbalanced=0.0001,
relStressTol=1.e-3,
doneHook='Finished()',
label='p3d'
),
# PyRunner(command='plotAddData()',iterPeriod=100),
]
def Finished():
O.save('post_iso.xml')
print 'Finished'
print (O.cell.size[0]*O.cell.size[1]*O.cell.size[2]-v)/v
print getStress()
print unbalancedForce()
O.pause()
O.run()
########Undrained triaxial test#########
from yade import pack,plot,qt,export
from math import fabs
import numpy as np
O.load('post_iso.xml')
O.trackEnergy=True
setContactFriction(0.5)
filename1='Packing_volume.txt'
a=np.loadtxt(filename1)
Ep_ini=O.energy['elastPotential']
Edamp_ini=O.energy['nonviscDamp']
Ekin_ini=utils.kineticEnergy()
work=0
law=Law2_ScGeom_FrictPhys_CundallStrack(traceEnergy=True)
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[law]),
NewtonIntegrator(damping=0.2),
GlobalStiffnessTimeStepper(),
PyRunner(command='inputwork()',iterPeriod=1),
PyRunner(command='plotAddData()',iterPeriod=10),
]
def inputwork():
global work
dwork=(0.5*(O.cell.prevVelGrad+O.cell.velGrad)*utils.getStress()).trace()*O.dt*(O.cell.hSize).determinant()
work=work+dwork
def plotAddData():
global Ep_ini
global Edamp_ini
global Ekin_ini
global work
plot.addData(
iter=O.iter,iter_=O.iter,
sxx=utils.getStress()[0,0],
syy=utils.getStress()[1,1],
szz=utils.getStress()[2,2],
exx=O.cell.getSmallStrain()[0,0],
eyy=O.cell.getSmallStrain()[1,1],
ezz=O.cell.getSmallStrain()[2,2],
exx1=O.cell.size[0],
eyy1=O.cell.size[1],
ezz1=O.cell.size[2],
ex1=O.cell.hSize[0,0],
ey1=O.cell.hSize[1,1],
ez1=O.cell.hSize[2,2],
exy1=O.cell.hSize[0,1],
eyx1=O.cell.hSize[1,0],
Z=avgNumInteractions(),
Zm=avgNumInteractions(skipFree=True),
voidratio=(O.cell.size[0]*O.cell.size[1]*O.cell.size[2]-a)/a,
# poros=voxelPorosity(500,(0,0,0),O.cell.size),
unbalanced=utils.unbalancedForce(),
t=O.time,
inwork=work
gWork=O.energy['gravWork'],
# Ep=O.energy['elastPotential'],
Ep=law.elasticEnergy()-Ep_ini,
Edamp=O.energy['nonviscDamp']-Edamp_ini,
Ediss=law.plasticDissipation(),
# Ediss=O.energy['plastDissip'],
Ekin=utils.kineticEnergy(),
Etot=law.elasticEnergy()-Ep_ini+O.energy['nonviscDamp']-Edamp_ini+law.plasticDissipation()+utils.kineticEnergy()-Ekin_ini,
# Etot=O.energy.total(),**O.energy
)
plot.saveDataTxt('macroFile',vars=('t','exx','eyy','ezz','sxx','syy','szz','Z','Zm','voidratio'))
plot.saveDataTxt('energyFile',vars=('t','unbalanced','gWork','Edamp','Ekin','Ep','Ediss','Etot','inwork'))
O.cell.trsf=Matrix3.Identity
O.cell.velGrad=Matrix3(-0.1,0,0,0,0.05,0,0,0,0.05)
O.run()
--
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
