yade-users team mailing list archive

[Horation <question703667@xxxxxxxxxxxxxxxxxxxxx>]

New question #703667 on Yade:
https://answers.launchpad.net/yade/+question/703667

I want to conduct an anisotropic consolidation test. The condition is that sigma3/sigma1=0.2 keeps constant, and during the consolidation process, the sigma1 and sigma3 increase linearly and their increments can be different but constant.Eventually sigma1 reaches 1000kPa and sigma3 reaches 200kPa. Here is my current code, I try to divide the process into many stages , but the data is not satisfying . Can anyone help me to improve the code or offer other better methods?
# unicode: UTF-8
# For 2D biaxial simulation
# 21/10/2016
# Yade version 2016.06a-24-0557faf~trusty

#########################################
### Defining parameters and variables ###
#########################################

#Material constants 

FrictionAngle = 0

Damp = 0.25

#Confining variables
ConfPress1 = -1.0e6
ConfPress3 = -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
O.load('isocompact200kpa2.yade.bz2')

########################################
#import necessary packages
from yade import pack,plot,os,timing

triax=O.engines[4]

#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 = 0.01,
    stressMask = 7,
    goal1 = -10000,
    goal2 = -10000,
    goal3 = -10000,
    
)

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(),Bo1_Box_Aabb()]),
    InteractionLoop(
     [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
     [Ip2_FrictMat_FrictMat_MindlinPhys()],
     [Law2_ScGeom_MindlinPhys_Mindlin()]
    ),
    GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
    triax,
    newton,
    PyRunner(iterPeriod=500,command='stressincrement()',label='dstress'),
    #PyRunner(realPeriod=2,command='endCheck()',label='check'),
    PyRunner(command='addPlotData()',iterPeriod=1000,label='record'),
    #VTKRecorder(iterPeriod=10000,recorders=['all'],fileName='./post/p1-'),
    TriaxialStateRecorder(iterPeriod=5000,file='WallStresses2')
]

mStress=(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3.0
#triax.max_vel=2*abs(-10000-mStress)/10000
# Simulation stop conditions defination
times = 0
def stressincrement():
    global times
    s1 = triax.stress(triax.wall_right_id)[0]
    s2 = triax.stress(triax.wall_top_id)[1]
    s3 = triax.stress(triax.wall_front_id)[2]
    unb=unbalancedForce()
    if abs(triax.goal1-s1)/(-triax.goal1)<0.01 and abs(triax.goal2-s2)/(-triax.goal2)<0.01 and abs(triax.goal3-s3)/(-triax.goal3)<0.01:
         times = times + 1
         triax.goal1 = triax.goal3 = -10000-times*100
         triax.goal2 = 2*triax.goal3
         print ('triax.goal1: %f, triax.goal2: %f, times: %f' %(triax.goal1,triax.goal2,times))
    if triax.goal1 < -200000:
         triax.goal1 = triax.goal3 = -200000
         triax.goal2 = -1000000
         if unb<0.05 and abs(ConfPress3-s1)/(-ConfPress3)<0.01 and abs(ConfPress1-s2)/(-ConfPress1)<0.01 and abs(ConfPress3-s3)/(-ConfPress3)<0.01:
             O.pause()
             O.save('isocompact200kpa2.yade.bz2')
#def endCheck():
    #unb=unbalancedForce()
    #s1 = triax.stress(triax.wall_right_id)[0]
    #s2 = triax.stress(triax.wall_top_id)[1]
    #s3 = triax.stress(triax.wall_front_id)[2]
    #if unb<0.05 and abs(ConfPress1-s1)/(-ConfPress1)<0.01 and abs(ConfPress3-s2)/(-ConfPress3)<0.01 and abs(ConfPress3-s3)/(-ConfPress3)<0.01:
       #O.pause()
       #O.save('isocompact200kpa2.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, disx: %f, disy: %f, disz: %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,triax.width,triax.height,triax.depth)) 
        #Output()
    ini_e22b = e22     
    plot.saveDataTxt('loadinglog2.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.