New question #706655 on Yade:
https://answers.launchpad.net/yade/+question/706655
I want to model a triaxial test on soft rock using Yade. I generate a sample
using Hertz-Mindlin law, and the sample contains coarse and fine grains
(specified PSD). Then in triaxial shearing stage, I change the material and
contact law into JCF material and contact law. But failed. This way is
impossible in Yade? Could you give me some advice? Thanks!
from __future__ import print_function
from yade import pack,plot,os,timing
from yade import export, ymport
import matplotlib;matplotlib.rc('axes',grid=True)
import pylab
import numpy as np
O.load('isocompact200kpa2.yade.bz2')
#O.reset()
for b in O.bodies:
if isinstance(b.shape, Sphere):
O.bodies.erase(b.id) #delete grains firstly and they will be
generate later
density=2600
frictionangle=0
poisson=0.15
young=30e9
damp=0.25
number=10000
walldensity=0
wallfrictionangle=0
wallpoisson=0.5
wallyoung=30000e9
mn=Vector3(0.,0.,0.)
mx=Vector3(2,2,2)
conpress=-200000
rate=0.1
#mat=O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=radians(frictionangle),density=density,label='spheres'))
wallmat=O.materials.append(JCFpmMat(type=1,frictionAngle=0,density=0,label='walls'))
idRockTest =
grainmat=O.materials.append(JCFpmMat(type=1,young=30e9,density=2500.0,poisson=0.1,frictionAngle=np.radians(18.0),tensileStrength=1e6,cohesion=1e6,label='Rock'))
O.bodies[0].material=O.materials[wallmat]
O.bodies[1].material=O.materials[wallmat]
O.bodies[2].material=O.materials[wallmat]
O.bodies[3].material=O.materials[wallmat]
O.bodies[4].material=O.materials[wallmat]
O.bodies[5].material=O.materials[wallmat]
#wallIds=O.bodies.append(aabbWalls([mn,mx],thickness=0.001,material='walls'))
#Material constants
sp3 = ymport.text('/home/hezihao/desk/bond/cloud.txt')
for s in range(0,len(sp3)):
#sp2[s].shape.color = (1.0,1.0,0)
sp3[s].state.material=grainmat
O.bodies.append(sp3[s])
for b in O.bodies:
if isinstance(b.shape,Sphere):
r=b.shape.radius
if r>0.03 and r<0.50:
b.shape.color=(1.0,0/255.,51/255.)
#b.mat=O.materials[grainmat]
if r>0.0 and r<0.03:
b.shape.color=(1.0,1.0,0)
#b.mat=O.materials[grainmat]
#for i in range(5,number+1):
#O.bodies[i].material=O.materials[grainmat]
#FrictionAngle = 35
Damp = 0.25
#Confining variables
ConfPress = -2.0e5
#Loading control
LoadRate = -0.01
#folders for post-processing
#try:
#os.mkdir('post')
#except:
#pass
#try:
#os.mkdir('pdata')
#except:
#pass
#try:
#os.mkdir('cdata')
#except:
#pass
#restart
########################################
#import necessary packages
from yade import pack,plot,os,timing
#def Output():
#e22=-triax.strain[1]
#particle info
#f=open('./pdata/pInfo'+'{:.5f}'.format(e22), 'w')
#f.write('Box position: \n')
#f.write('leftwall:%.5f %.5f
%.5f\n'%(O.bodies[0].state.pos[0],O.bodies[0].state.pos[1],O.bodies[0].state.pos[2]))
#f.write('rightwall:%.5f %.5f
%.5f\n'%(O.bodies[1].state.pos[0],O.bodies[1].state.pos[1],O.bodies[1].state.pos[2]))
#f.write('topwall:%.5f %.5f
%.5f\n'%(O.bodies[3].state.pos[0],O.bodies[3].state.pos[1],O.bodies[3].state.pos[2]))
#f.write('bottomwall:%.5f %.5f
%.5f\n'%(O.bodies[2].state.pos[0],O.bodies[2].state.pos[1],O.bodies[2].state.pos[2]))
#f.write('frontwall:%.5f %.5f
%.5f\n'%(O.bodies[5].state.pos[0],O.bodies[5].state.pos[1],O.bodies[5].state.pos[2]))
#f.write('backwall:%.5f %.5f
%.5f\n'%(O.bodies[4].state.pos[0],O.bodies[4].state.pos[1],O.bodies[4].state.pos[2]))
#f.write('ID x y z radius disx disy disz rotx roty rotz \n')
#for b in O.bodies:
#if isinstance(b.shape, Sphere):
#pos=b.state.pos
#rad=b.shape.radius
#displ=b.state.displ()
#rot=b.state.rot()
#f.write('%-16g %-16g %-16g %-16g %-16g %-16g %-16g %-16g %-16g %-16g
%-16g\n'%(b.id,pos[0],pos[1],pos[2],rad,displ[0],displ[1],displ[2],rot[0],rot[1],rot[2]))
#f.close()
#contact info
#g = open('./cdata/cInfo'+'{:.5f}'.format(e22), 'w')
#print ('Contact information at Iter %d' % O.iter)
#g.write('ctype id1 id2 nfx nfy nfz tfx tfy tfz \n')
#for k in O.interactions:
#if isinstance(O.bodies[k.id1].shape,Sphere) and
isinstance(O.bodies[k.id2].shape,Box):
#g.write('01:%.5f %.5f %.5f %.5f %.5f %.5f %.5f
%.5f\n'%(k.id1,k.id2,k.phys.normalForce[0],
k.phys.normalForce[1],k.phys.normalForce[2],k.phys.shearForce[0],k.phys.shearForce[1],k.phys.shearForce[2]))
#elif isinstance(O.bodies[k.id1].shape,Box) and
isinstance(O.bodies[k.id2].shape,Sphere):
#g.write ('10:%.5f %.5f %.5f %.5f %.5f %.5f %.5f
%.5f\n'%(k.id1,k.id2,k.phys.normalForce[0],
k.phys.normalForce[1],k.phys.normalForce[2],k.phys.shearForce[0],k.phys.shearForce[1],k.phys.shearForce[2]))
#elif isinstance(O.bodies[k.id1].shape,Sphere) and
isinstance(O.bodies[k.id2].shape,Sphere):
#g.write ('00:%.5f %.5f %.5f %.5f %.5f %.5f %.5f
%.5f\n'%(k.id1,k.id2,k.phys.normalForce[0],
k.phys.normalForce[1],k.phys.normalForce[2],k.phys.shearForce[0],k.phys.shearForce[1],k.phys.shearForce[2]))
#g.close()
#Output()
#################################
#####Defining triaxil engines####
#################################
triax=TriaxialStressController(
#define wall ids
#wall_bottom_id = wallIds[2],
#wall_top_id = wallIds[3],
#wall_left_id = wallIds[0],
#wall_right_id = wallIds[1],
#wall_back_id = wallIds[4],
#wall_front_id = wallIds[5],
thickness = 0.001,
#maxMultiplier = 1.002,
internalCompaction = False, # If true the confining pressure is generated
by growing particles
max_vel = 1.0,
stressMask = 5,
goal1 = ConfPress,
goal2 = LoadRate,
goal3 = ConfPress,
)
ini_e22a = -triax.strain[1]
ini_e22b = -triax.strain[1]
O.usesTimeStepper=True
#O.dt=2e-6
newton=NewtonIntegrator(damping=Damp)
#setContactFriction(radians(FrictionAngle))
O.trackEnergy=True
O.timingEnabled=True
###engine
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.3),Bo1_Box_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.3),Ig2_Box_Sphere_ScGeom()],
[Ip2_JCFpmMat_JCFpmMat_JCFpmPhys()],
[Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(label='interactionLaw')]
),
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
triax,
newton,
PyRunner(realPeriod=2,command='endCheck()',label='check'),
PyRunner(command='addPlotData()',iterPeriod=100,label='record'),
#VTKRecorder(iterPeriod=100,recorders=['all'],fileName='./post/p1-'),
#TriaxialStateRecorder(iterPeriod=10000,file='WallStresses3')
]
# Simulation stop conditions defination
def endCheck():
unb=unbalancedForce()
e22=-triax.strain[1]
if abs(e22-0.40)<0.001:
O.pause()
O.save('e0.05.yade.bz2')
# collect history of data
def addPlotData():
global ini_e22a
global ini_e22b
e22=-triax.strain[1]
unb = unbalancedForce()
mStress =
-(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3.0
volume =
(O.bodies[1].state.pos[0]-O.bodies[0].state.pos[0])*(O.bodies[3].state.pos[1]-O.bodies[2].state.pos[1])*(O.bodies[5].state.pos[2]-O.bodies[4].state.pos[2])
Porosity=1-sum(4*pi*b.shape.radius**3/3 for b in O.bodies if
isinstance(b.shape,Sphere))/volume
plot.addData(e11=-triax.strain[0], e22=-triax.strain[1],
e33=-triax.strain[2],
ev=-triax.strain[0]-triax.strain[1]-triax.strain[2],
s11=-triax.stress(triax.wall_right_id)[0],
s22=-triax.stress(triax.wall_top_id)[1],
s33=-triax.stress(triax.wall_front_id)[2],
ub=unbalancedForce(),
dstress=(-triax.stress(triax.wall_top_id)[1])-(-triax.stress(triax.wall_right_id)[0]),
disx=triax.width,
disy=triax.height,
disz=triax.depth,
i=O.iter,
porosity=Porosity,
cnumber=avgNumInteractions(),
ke=utils.kineticEnergy(),
totale=O.energy.total()
)
#if (e22 - ini_e22a) > 0.001 and e22 < 0.05:
# O.save('save'+'{:.4f}'.format(e22)+'yade.bz2')
# ini_e22a = e22
if (e22 - ini_e22b) > 0.00005:
print ('unbalanced force: %f, mean stress: %f, s11: %f, s22: %f, s33:
%f, coordination number: %f, porosity: %f'
%(unb,mStress,-triax.stress(triax.wall_right_id)[0],-triax.stress(triax.wall_top_id)[1],-triax.stress(triax.wall_front_id)[2],avgNumInteractions(),Porosity))
#Output()
ini_e22b = e22
plot.saveDataTxt('loadinglog3.txt.bz2')
#plot.saveGnuplot('plotScript')
##################################################
######Defining output for postprocessing##########
##################################################
###display
#yade.qt.Controller(),yade.qt.View()
### declare what is to plot. "None" is for separating y and y2 axis
#plot.plots={'i':('ub',),'i ':('e11','e22','e33',),'
i':('s11','s22','s33',),'e22':'dstress',' e22':'ev'}
#plot.plots={'i':('ub',),'i ':('s11','s22','s33'),' i':('e11','e22','e33')}
### the traditional triaxial curves would be more like this:
plot.plots={'e22':('s22',None,'ev')}
plot.plot()
#O.run()#; O.wait()
--
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
