Question #695935 on Yade changed:
https://answers.launchpad.net/yade/+question/695935
Yuxuan Wen gave more information on the question:
I run this file first:
#-----------------------------------------------------------------------------------------------------
## length (m), time (s), mass (1kg), force (N), pressure (Pa)
from yade import pack, plot, qt, export, os
##
lx=0.2
ly=0.2
lz=0.14
target_r=0.005 # target r=0.005m=0.5cm, d=0.01m=1cm, l=20d
target_d=2*target_r
target_phi=0.550
target_e=(1-target_phi)/target_phi
num_spheres=int(lx*ly*lz*target_phi/(4.0/3.0*3.1415926*target_r*target_r*target_r))
flag=1 # mark time when changing friction
## space for generating particle cloud
thickness=target_d # top and bottom walls' thickness
wallsurfh=target_r # distance between wall surface to the cell
mn=Vector3(target_d,target_d,wallsurfh+thickness+target_d)
mx=Vector3(2*lx-target_d,2*ly-target_d,wallsurfh+thickness+2*lz-target_d)
## material properties
density=2650.0 # kg/m3
kn=4e5 # kn/2 is the stiffness of the contact
ks=kn*2.0/7.0
gamma=50.0
cn=4.0/3.0*3.1415926*(0.005*0.005*0.005)*density/4*gamma*2 #0.03469
frict=0
frictnew=26.57 # after stable, change the particle from frictionless to
frictional
## insert periodic cell
O.periodic=True
O.cell.hSize=Matrix3( 2*lx, 0, 0,
0, 2*ly, 0,
0, 0, wallsurfh+thickness+2*lz+thickness+wallsurfh)
## insert the walls for consolidation
O.materials.append(ViscElMat(kn=kn,ks=ks,frictionAngle=radians(frict),cn=cn,cs=0,density=density,label='Box'))
lowbox=O.bodies.append(utils.box(center=(lx,ly,wallsurfh+0.5*thickness),
extents=(5,5,0.5*thickness), fixed=False, material='Box',wire=False))
upbox=O.bodies.append(utils.box(center=(lx,ly,wallsurfh+thickness+2*lz+0.5*thickness),
extents=(5,5,0.5*thickness), fixed=False, material='Box',wire=False))
## use a SpherePack object to generate a random loose particles packing
O.materials.append(ViscElMat(kn=kn,ks=ks,frictionAngle=radians(frict),cn=cn,cs=0,density=density,label='spheres'))
sp=pack.SpherePack()
sp.makeCloud(mn,mx,rMean=target_r,num=num_spheres,periodic=True,seed=1) #"seed"
make the "random" generation always the same
O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])
msphere=0.0
sphereid=[]
for i in O.bodies:
if i.id>=2:
msphere=msphere+i.state.mass
sphereid.append(i.id)
######################### define engines #########################
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()],allowBiggerThanPeriod=True),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
[Law2_ScGeom_ViscElPhys_Basic()]
),
PyRunner(command='servo()',iterPeriod=1), # add force must before
NewtonIntegrator(), otherwise it will be ForceResetter()
GlobalStiffnessTimeStepper(),
NewtonIntegrator(damping=0.2), # gravity=(0,0,0), numerical damping,
dissipate Ekine because of quasi-static loading
PyRunner(command='addPlotData()',iterPeriod=200),
PyRunner(command='finished()',iterPeriod=200),
#PyRunner(command='fabricdata()',iterPeriod=500),
#VTKRecorder(fileName='vtkfile/consolidation-',recorders=['spheres','intr','velocity','stress','force'],iterPeriod=2000),
#qt.SnapshotEngine(iterPeriod=2000,fileBase='video/consolidation-',label='snapshooter')
]
#plot.live=True
plot.plots={'t':('sxx','syy','szz'),'t ':('exx','eyy','ezz'),'t ':('void'),'t
':('Eunbal')}
plot.plot()
#O.run()
######################### define functions #########################
def servo():
O.bodies[upbox].dynamic=False
O.bodies[lowbox].dynamic=False
if O.cell.size[0] > (lx+1e-8):
rate1=0.1
else:
rate1=0
O.cell.velGrad=Matrix3(-rate1,0,0, 0,-rate1,0, 0,0,0)
if O.bodies[upbox].state.pos[2]-O.bodies[lowbox].state.pos[2]-thickness
> (lz+1e-8):
rate2=0.1
else:
rate2=0
vz=rate2*(O.bodies[upbox].state.pos[2]-O.bodies[lowbox].state.pos[2]-thickness)/2
O.bodies[lowbox].state.vel=(0,0,vz)
O.bodies[upbox].state.vel=(0,0,-vz)
def addPlotData():
V=O.cell.size[0]*O.cell.size[1]*(O.bodies[upbox].state.pos[2]-O.bodies[lowbox].state.pos[2]-thickness)
Vs=msphere/density
n0=0
n1=0
for i in sphereid:
if len(O.bodies[i].intrs())==0:
n0=n0+1
if len(O.bodies[i].intrs())==1:
n1=n1+1
plot.addData(
nsphere=num_spheres,
N0=n0,
N1=n1,
step=O.iter,
t=O.time,
phi=Vs/V,
void=(V-Vs)/Vs,
sxx=utils.getStress(volume=V)[0][0],
sxy=utils.getStress(volume=V)[0][1], sxz=utils.getStress(volume=V)[0][2],
syx=utils.getStress(volume=V)[1][0],
syy=utils.getStress(volume=V)[1][1], syz=utils.getStress(volume=V)[1][2],
szx=utils.getStress(volume=V)[2][0],
szy=utils.getStress(volume=V)[2][1], szz=utils.getStress(volume=V)[2][2],
exx=O.cell.size[0],
eyy=O.cell.size[1],
ezz=O.bodies[upbox].state.pos[2]-O.bodies[lowbox].state.pos[2]-thickness,
s_up=O.forces.f(upbox)[2]/(O.cell.size[0]*O.cell.size[1]),
s_low=O.forces.f(lowbox)[2]/(O.cell.size[0]*O.cell.size[1]),
Z=avgNumInteractions(), # interactions between clump
particle-walls are also considered in YADE 2018.02b
Zm=avgNumInteractions(skipFree=True),
Ekine=utils.kineticEnergy(),
Etot=utils.O.energy.total,
Eunbal=utils.unbalancedForce(),**O.energy
# **O.energy converts dictionary-like O.energy to plot.addData
arguments
)
plot.saveDataTxt(os.getcwd()+'/normalfiles/consolidation',vars=('t','step','sxx','syy','szz','exx','eyy','ezz','phi','nsphere','Eunbal','N0','N1'))
# overall state variables under different step are saved in one file, no loop
def fabricdata():
filename=os.getcwd()+'/fabricdata/'+'consolidation'+str(O.iter)
f = open(filename,'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
point = i.geom.contactPoint
norm = i.geom.normal
ovp = i.geom.penetrationDepth
Fn = i.phys.normalForce
Fs = i.phys.shearForce
f.write('%-8i %-8i %-12g %-12g %-12g %-12g %-12g %-12g %-12g
%-12g %-12g %-12g %-12g %-12g
%-12g\n'%(i.id1,i.id2,point[0],point[1],point[2],norm[0],norm[1],norm[2],Fn[0],Fn[1],Fn[2],Fs[0],Fs[1],Fs[2],ovp))
f.close
def finished():
V=O.cell.size[0]*O.cell.size[1]*(O.bodies[upbox].state.pos[2]-O.bodies[lowbox].state.pos[2]-thickness)
Vs=msphere/density
phi=Vs/V
void=(V-Vs)/Vs
Eunbal=utils.unbalancedForce()
global flag, step1
if (Eunbal<=1e-6 or
O.time>=1e4*sqrt(4.0/3.0*3.1415926*target_r*target_r*target_r*density/(kn/2)))
and flag==1:
flag=2
step1=O.iter
setContactFriction(radians(frictnew)) # frictionAngle of boxes
are still 0 as they are not dynamic body
print ('friction added at '+ str(O.iter))
print ('step='+str(O.iter))
print ('phi='+str(phi),
'error='+str(fabs((target_phi-phi)/target_phi)))
print ('void ratio='+str(void),
'error='+str(fabs((target_e-void)/target_e)))
print ('radius='+str(O.bodies[2].shape.radius))
print ('Unbalanced force = '+str(Eunbal))
elif (Eunbal<=1e-6 or
O.time>=1e4*sqrt(4.0/3.0*3.1415926*target_r*target_r*target_r*density/(kn/2)))
and O.iter-step1>=10000:
fabricdata()
O.save('consolidation.xml')
print ('Consolidation Finished')
print ('step='+str(O.iter))
print ('phi='+str(phi),
'error='+str(fabs((target_phi-phi)/target_phi)))
print ('void ratio='+str(void),
'error='+str(fabs((target_e-void)/target_e)))
print ('radius='+str(O.bodies[2].shape.radius))
print ('Unbalanced force = '+str(Eunbal))
O.pause()
--
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 : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
