Hiya All, I wrote a code that is supposed to simulate a Triaxial_Undrained test of a saturated sphere packing within a cylindrical specimen. The The code run without any issue but it's not working properly. Can anyone help me to make this code work? The code is copied below: # import corrosponding yade modulus from yade import pack from yade import geom from yade import plot from yade import ymport from yade import qt from yade.gridpfacet import * import gts, os.path, locale locale.setlocale(locale.LC_ALL, 'en_US.UTF-8') import sys,os from yade import plot,pack,export sys.path.append(os.pardir) ######################################################################## # A. DEFINING VARIABLES, MATERIALS' PROPERTIES AND PACKING # ######################################################################## # A.a). define variables key = 'Triaxial_Undrained' # file name to be saved young=550e6 # normal contact stiffness compFricDegree = 1.8 # initial contact friction during the confining phase finalFricDegree = 43 # contact friction during the deviatoric loading poisson = 0.3 # shear-to-normal stiffness ratio isoStress = 110000 # confining stress conStress = 100000 # confining stress for deviatoric loading stage width = 1.4e-1 # sample width height = 2.8e-1 # target sample height(after consolidation) height_0 = 3.2e-1 # initial sample height num_spheres=500 # number of spheres R_p = 0.0084 # mean particle radius rCoff = 10 # thickness of top and bot sphere cap (based on rParticle) rParticle = 0.02e-1 # membrane grid seed size alpha = 8 rate = 0.1 # loading rate (strain rate) damp = 0.3 # damping coefficient targetPorosity = 0.43 # target porosity thresholdvalue = 0.05 # threshold unbalance force final_rate = 0.1 # strain rate for deviator loading thresholdstrain = 0.06 # threshold axial strain for terminate enlargefactor = 1.00 # A.b). create materials for sand spheres and plates Sand = O.materials.append(CohFrictMat( young=young,poisson=poisson,frictionAngle=radians(compFricDegree), alphaKr=0.25,alphaKtw=0,etaRoll=0.005,etaTwist=0, normalCohesion=5e6,shearCohesion=5e6, momentRotationLaw=True,density=2650,label='spheres' )) # A.c). create membrane materials 207 GridMat = O.materials.append(CohFrictMat( young=100e6,poisson=0.3,density=2650,frictionAngle=radians(0), alphaKr=0,alphaKtw=0,etaRoll=0,etaTwist=0, normalCohesion=1e9,shearCohesion=1e9, momentRotationLaw=True,label='gridNodeMat' )) pFacetMat = O.materials.append(FrictMat( young=100e6,poisson=0.3,density=2650,frictionAngle=radians(0),label='pFacetMat' )) # A.d). create TOP & BOT plate materials frictMat = O.materials.append(FrictMat( young=100e6,poisson=0.3,density=2650,frictionAngle=radians(0),label='frictMat' )) ################################################################### # B. DEFINING GLOBAL ENGINES # ################################################################### #**********************************************************************# O.engines=[ ForceResetter(), InsertionSortCollider([ Bo1_Sphere_Aabb(), Bo1_PFacet_Aabb(), Bo1_Facet_Aabb(), Bo1_GridConnection_Aabb() ]), InteractionLoop( [ Ig2_Sphere_Sphere_ScGeom6D(), Ig2_GridNode_GridNode_GridNodeGeom6D(), Ig2_GridConnection_GridConnection_GridCoGridCoGeom(), Ig2_Sphere_PFacet_ScGridCoGeom(), Ig2_Facet_Sphere_ScGeom6D() ], [ Ip2_FrictMat_FrictMat_FrictPhys(), Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label="cohesiveIp") ], [ Law2_ScGeom_FrictPhys_CundallStrack(), Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm =True,always_use_moment_law=False,label='cohesiveLaw'), Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(), Law2_ScGridCoGeom_FrictPhys_CundallStrack(), ],label="iloop" ), GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=25,timestepSafetyCoefficient=0.8), NewtonIntegrator(gravity=(0,0,0),damping=0.1,label='newton') ] #**********************************************************************# O.engines[2].lawDispatcher.functors[1].always_use_moment_law=False O.engines[2].physDispatcher.functors[1].setCohesionOnNewContacts=False ####################################################################### # C. GENERATING PACKING # ####################################################################### # C.a). generate random dense sphere pack sp=pack.SpherePack() pred = pack.inCylinder((0,0,0),(0,0,height_0),.5*width) sp.makeCloud((0,0,0),(.3,.3,2),rMean=0.0083,rRelFuzz=0.1) #sp = pack.randomDensePack( pred,spheresInCell=num_spheres,radius=R_p,rRelFuzz=0.3, returnSpherePack=True,memoDbg=True),#memoizeDb='/tmp/loosePackings11.sqlite') sp.toSimulation(color=(0,1,1),material=Sand) triax=TriaxialStressController( maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth) finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth) thickness = 0, stressMask = 7, max_vel = 0.005, internalCompaction=True, # If true the confining pressure is generated by growing particles ) newton=NewtonIntegrator(damping=0.2) O.engines=[ ForceResetter(), InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]), InteractionLoop( [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()], [Ip2_FrictMat_FrictMat_FrictPhys()], [Law2_ScGeom_FrictPhys_CundallStrack()],label="iloop" ), FlowEngine(dead=1,label="flow"),#introduced as a dead engine for the moment, see 2nd section GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8), triax, newton ] triax.goal1=triax.goal2=triax.goal3=-10000 while 1: O.run(1000, True) unb=unbalancedForce() if unb<0.001 and abs(-10000-triax.meanStress)/10000<0.001: break setContactFriction(radians(finalFricDegree)) sand=sp.toSimulation(color=(0,1,1),material=Sand) # C.b). define different sections of sphere pack bot = [O.bodies[s] for s in sand if O.bodies[s].state.pos[2]height_0-rParticle*rCoff] tot = [O.bodies[s] for s in sand if O.bodies[s].state.pos[2]<=height_0] # C.c). detect the position of particles in top & bot layer top_limit = 0 top_id = 0 for s in top: if s.state.pos[2]>=top_limit: top_limit = s.state.pos[2] top_id = s.id bot_limit = height_0 bot_id = 0 for s in bot: if s.state.pos[2]<=bot_limit: bot_limit = s.state.pos[2] bot_id = s.id # C.d). create facet wall around particle packing facets = [] nw = 45 nh = 15 rCyl2 = .5*width / cos(pi/float(nw)) for r in xrange(nw): for h in xrange(nh): v1 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw)), rCyl2*sin(2*pi*(r+0)/float(nw)), height_0*(h+0)/float(nh) ) v2 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw)), rCyl2*sin(2*pi*(r+1)/float(nw)), height_0*(h+0)/float(nh) ) v3 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw)), rCyl2*sin(2*pi*(r+1)/float(nw)), height_0*(h+1)/float(nh) ) v4 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw)), rCyl2*sin(2*pi*(r+0)/float(nw)), height_0*(h+1)/float(nh) ) f1 = facet((v1,v2,v3),color=(0,0,1),material=frictMat) f2 = facet((v1,v3,v4),color=(0,0,1),material=frictMat) facets.extend((f1,f2)) wall = O.bodies.append(facets) # C.e). define different sections of facet wall for b in wall: O.bodies[b].state.blockedDOFs = 'xyzXYZ' O.bodies[b].state.vel = (0,0,0) # C.f). create top & bot facet plate facets3 = [] nw=45 rCyl2 = (.6*width+2*rParticle) / cos(pi/float(nw)) for r in xrange(nw): if r%2==0: v1 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw)), rCyl2*sin(2*pi*(r+0)/float(nw)), height_0 ) v2 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw)), rCyl2*sin(2*pi*(r+1)/float(nw)), height_0 ) v3 = Vector3( rCyl2*cos(2*pi*(r+2)/float(nw)), rCyl2*sin(2*pi*(r+2)/float(nw)), height_0 ) v4 = Vector3( 0, 0, height_0 ) f1 = facet((v1,v2,v4),color=(0,0,0),material=frictMat) f2 = facet((v2,v3,v4),color=(0,0,0),material=frictMat) facets3.extend((f1,f2)) topcap = O.bodies.append(facets3) facets3 = [] for r in xrange(nw): if r%2==0: v1 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw)), rCyl2*sin(2*pi*(r+0)/float(nw)), 0 ) v2 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw)), rCyl2*sin(2*pi*(r+1)/float(nw)), 0 ) v3 = Vector3( rCyl2*cos(2*pi*(r+2)/float(nw)), rCyl2*sin(2*pi*(r+2)/float(nw)), 0 ) v4 = Vector3( 0, 0, 0 ) f1 = facet((v1,v2,v4),color=(0,0,0),material=frictMat) f2 = facet((v2,v3,v4),color=(0,0,0),material=frictMat) facets3.extend((f1,f2)) botcap = O.bodies.append(facets3) # C.g). define top & bot wall id for s in topcap: top_id = s bot_id = 0 for s in botcap: bot_id = s # D.h). calculate porosity V_sand = 0 num_sand = 0 for b in sand: r = O.bodies[b].shape.radius V_sand += 1.3333333*3.1416*r*r*r num_sand +=1 porosity = 1-V_sand/(.25*width*width*3.1416*height_0) print 'v_sand= ',V_sand,' number of sand: ',num_sand,'porosity is: ',porosity O.pause() triax.stressMask=2 triax.goal1=triax.goal3=0 triax.internalCompaction=False triax.wall_bottom_activated=False #load triax.goal2=-11000; O.run(2000,1) #unload triax.goal2=-10000; O.run(2000,1) #load triax.goal2=-11000; O.run(2000,1) e22=triax.strain[1] #unload triax.goal2=-10000; O.run(2000,1) e22=e22-triax.strain[1] modulus = 1000./abs(e22) ####################################################################### # D. DEFINING ADD-ON FUNCTIONS # ####################################################################### # D.a). a function for saving variables def plotAddData(): f1 = sum(O.forces.f(b)[2] for b in topcap) f2 = sum(O.forces.f(b)[2] for b in botcap) f11 = sum(O.forces.f(b.id)[2] for b in top) f22 = sum(O.forces.f(b.id)[2] for b in bot) fa = abs(.5*(f2-f1)) fa1 = abs(.5*(f22-f11)) e = (top[0].state.displ()[2] - bot[0].state.displ()[2])/height_0-e_ini f4 = 0 r_cum = 0 count = 0 Area = 0 for f in membrane_grid: f.shape.color = Vector3(0,0,0) x = f.state.pos[0] y = f.state.pos[1] dist = math.sqrt(x*x+y*y) n = Vector3(x/dist,y/dist,0) z = f.state.pos[2] a = 2*3.1416*dist/(360/alpha)*6*rParticle*(1+e) if z<=O.bodies[top_id].state.pos[2] and z>=O.bodies[bot_id].state.pos[2]: count += 1 r_local = dist r_cum += r_local r_avg = r_cum/count-rParticle Area = r_avg*2*3.1416*(O.bodies[top_id].state.pos[2] - O.bodies[bot_id].state.pos[2]) area_avg = Area/count s = fa/(3.1416*r_avg*r_avg) s1 = fa1/(3.1416*r_avg*r_avg) for b in membrane_grid: x = b.state.pos[0] y = b.state.pos[1] dist = math.sqrt(x*x+y*y) n = Vector3(x/dist,y/dist,0) a = 2*3.1416*dist/(360/alpha)*6*rParticle*(1+e) z = b.state.pos[2] if z<=O.bodies[top_id].state.pos[2] and z>=O.bodies[bot_id].state.pos[2]: f_local = O.forces.f(b.id) length = math.sqrt(f_local[0]*f_local[0]+f_local[1]*f_local[1]+f_local[2]*f_local[2]) cos_theta = (n[0]*f_local[0]+n[1]*f_local[1]+n[2]*f_local[2])/length p_normal = (length*cos_theta/a) f4 += (p_normal) p = abs(f4/count/1000) h = O.bodies[top_id].state.pos[2] - O.bodies[bot_id].state.pos[2] VV = h*r_avg*r_avg*3.1416 dV = VV-V_ini ev = -((O.bodies[top_id].state.pos[2] - O.bodies[bot_id].state.pos[2])*r_avg*r_avg*3.1416-V_ini)/V_ini er = (r_avg-R_avg)/R_avg if (abs(e*100)>thresholdstrain*100): O.pause() plot.addData( i = O.iter, q = (abs(s)-p*1000)/1000, q1 = (abs(s1)-conStress)/1000, p = p, u = conStress/1000-p, e = -e*100, ev = ev*100, ) for b in tot: b.shape.color=scalarOnColorScale(b.state.rot().norm(),0,pi/2.) #return (dV,e) # D.b). a function for adding force (servo-controlled of lateral wall) def addforce(): h_sample = O.bodies[top_id].state.pos[2] - O.bodies[bot_id].state.pos[2] #print 'height is ',h_sample r_cum = 0 count = 0 f4 = 0 for b in membrane_grid: x = b.state.pos[0] y = b.state.pos[1] dist = math.sqrt(x*x+y*y) n = Vector3(x/dist,y/dist,0) a = 2*3.1416*dist/(360/alpha)*6*rParticle*(1+e) z = b.state.pos[2] if z<=O.bodies[top_id].state.pos[2] and z>=O.bodies[bot_id].state.pos[2]: f_local = O.forces.f(b.id) length = math.sqrt(f_local[0]*f_local[0]+f_local[1]*f_local[1]+f_local[2]*f_local[2]) cos_theta = (n[0]*f_local[0]+n[1]*f_local[1]+n[2]*f_local[2])/length p_normal = (length*cos_theta/a) f4 += (p_normal) count += 1 p = abs(f4/count/1000) for f in membrane_grid: f.shape.color = Vector3(0,0,0) x = f.state.pos[0] y = f.state.pos[1] dist = math.sqrt(x*x+y*y) n = Vector3(x/dist,y/dist,0) z = f.state.pos[2] a = 2*3.1416*dist/(360/alpha)*6*rParticle*(1+e) if z<=O.bodies[top_id].state.pos[2] and z>=O.bodies[bot_id].state.pos[2]: r_local = dist r_cum += r_local r_avg = r_cum/count-rParticle Volume = r_avg*r_avg*3.1416*h_sample delV = Volume - V_ini for b in topcap: O.bodies[b].state.blockedDOFs = 'xyzXYZ' O.bodies[b].state.vel = (0,0,-vel_a) for b in botcap: O.bodies[b].state.blockedDOFs = 'xyzXYZ' O.bodies[b].state.vel = (0,0,vel_a) for f in membrane_grid: f.shape.color = Vector3(0,0,0) x = f.state.pos[0] y = f.state.pos[1] dist = math.sqrt(x*x+y*y) n = Vector3(x/dist,y/dist,0) z = f.state.pos[2] a = 2*3.1416*dist/(360/alpha)*6*rParticle#*(1+e) if z<=O.bodies[top_id].state.pos[2] and z>=O.bodies[bot_id].state.pos[2]: f.state.blocked = 'z' f.state.vel = (dist/h_sample)*vel_a*n if delV>0: f.state.vel = -3.0*(dist/h_sample)*vel_a*n else: f.state.vel = 3.0*(dist/h_sample)*vel_a*n else: f.state.vel = 0*n # D.c). a function for recording data def checkrecord(): plot.saveDataTxt('results_'+key) # D.d). a function used for consolidation def confining(): e_ini = (top[0].state.displ()[2] - bot[0].state.displ()[2])/height_0 f1 = sum(O.forces.f(b)[2] for b in topcap) f2 = sum(O.forces.f(b)[2] for b in botcap) f4 = 0 r_cum = 0 count = 0 a = 2*3.1416*(.5*width+rParticle)/(360/alpha)*6*rParticle for f in membrane_grid: f.shape.color = Vector3(0,0,0) x = f.state.pos[0] y = f.state.pos[1] dist = math.sqrt(x*x+y*y) n = Vector3(x/dist,y/dist,0) a = 2*3.1416*dist/(360/alpha)*6*rParticle z = f.state.pos[2] f.state.vel = -vel_ini_r*n for b in topcap: O.bodies[b].state.blockedDOFs = 'xyzXYZ' O.bodies[b].state.vel = (0,0,-vel_ini_a) for b in botcap: O.bodies[b].state.blockedDOFs = 'xyzXYZ' O.bodies[b].state.vel = (0,0,vel_ini_a) for b in membrane_grid: x = b.state.pos[0] y = b.state.pos[1] dist = math.sqrt(x*x+y*y) n = Vector3(x/dist,y/dist,0) a = 2*3.1416*dist/(360/alpha)*6*rParticle z = b.state.pos[2] if z<=O.bodies[top_id].state.pos[2] and z>=O.bodies[bot_id].state.pos[2]: f_local = O.forces.f(b.id) length = math.sqrt(f_local[0]*f_local[0]+f_local[1]*f_local[1]+f_local[2]*f_local[2]) cos_theta = (n[0]*f_local[0]+n[1]*f_local[1]+n[2]*f_local[2])/length p_normal = (length*cos_theta/a) f4 += (p_normal) r_cum += dist count += 1 r_avg = r_cum/count-rParticle fa = abs(.5*(f2-f1)) p_r = f4/count/1000 p_a = fa/(3.1416*0.25*width*width)/1000 e_ini2 = ((O.bodies[top_id].state.pos[2] - O.bodies[bot_id].state.pos[2])- height_0)/height_0 V_ini = (O.bodies[top_id].state.pos[2] - O.bodies[bot_id].state.pos[2])*r_avg*r_avg*3.1416 R_avg = r_avg H_ini = O.bodies[top_id].state.pos[2] - O.bodies[bot_id].state.pos[2] porosity = 1-V_sand/((R_avg)*(R_avg)*3.1416*H_ini) #return (p_r,p_a,e_ini,V_ini,R_avg,H_ini,porosity) flow.dead=0 flow.defTolerance=0.3 flow.meshUpdateInterval=200 flow.useSolver=3 flow.permeabilityFactor=1 flow.viscosity=10 flow.bndCondIsPressure=[0,0,1,1,0,0] flow.bndCondValue=[0,0,1,0,0,0] flow.boundaryUseMaxMin=[0,0,0,0,0,0] O.dt=0.1e-3 O.dynDt=False O.run(1,1) Qin = flow.getBoundaryFlux(2) Qout = flow.getBoundaryFlux(3) permeability = abs(Qin)/1.e-4 #size is one, we compute K=V/H print "Qin=",Qin," Qout=",Qout," permeability=",permeability flow.bndCondIsPressure=[0,0,0,1,0,0] flow.bndCondValue=[0,0,0,0,0,0] newton.damping=0 Cv=permeability*modulus/1e4 zeroTime=O.time zeroe22 = - triax.strain[1] dryFraction=0.05 #the top layer is affected by drainage on a certain depth, we account for it here drye22 = 1000/modulus*dryFraction wetHeight=1*(1-dryFraction) # D.e). a function for stablization def stable(): for f in membrane_grid: x = f.state.pos[0] y = f.state.pos[1] dist = math.sqrt(x*x+y*y) n = Vector3(x/dist,y/dist,0) a = 2*3.1416*dist/(360/alpha)*6*rParticle z = f.state.pos[2] f.state.blockedDOFs = 'xyzXYZ' f.state.vel = 0*n for b in topcap: O.bodies[b].state.blockedDOFs = 'xyzXYZ' O.bodies[b].state.vel = (0,0,0) for b in botcap: O.bodies[b].state.blockedDOFs = 'xyzXYZ' O.bodies[b].state.vel = (0,0,0) def compress(): for b in wall: O.bodies[b].state.blockedDOFs = 'xyzXYZ' O.bodies[b].state.vel = (0,0,0) ####################################################################### # E. APPLYING CONFINING STRESS TO FLEXIBLE MEMBRANE # ####################################################################### # E.a). adding corrosponding python function O.engines=O.engines+[ PyRunner(command='compress()',iterPeriod=1), PyRunner(command='plotAddData',iterPeriod=1) ] # E.b). compress until target porosity vel_ini_a = rate*height_0 vel_ini_r = rate*height_0 for b in topcap: O.bodies[b].state.blockedDOFs = 'xyzXYZ' O.bodies[b].state.vel = (0,0,-vel_ini_a) for b in botcap: O.bodies[b].state.blockedDOFs = 'xyzXYZ' O.bodies[b].state.vel = (0,0,vel_ini_a) while 1: O.run(100,True) h = (O.bodies[top_id].state.pos[2]-O.bodies[bot_id].state.pos[2]) V = h*0.25*width*width*3.1416 porosity = 1-V_sand/V print 'porosity: ',porosity, ' height: ', h if (porosity <= targetPorosity): print 'compression stage finished!' break h_ini = h # record height after compression # E.c). ADD CONFINING STRESS # E.c.1). remove facet wall for b in wall: O.bodies.erase(b) # E.c.2). create membrane around particle packing by reading mesh file shiftfactor = O.bodies[bot_id].state.pos[2]-((height-h_ini)/2) nodesIds,cylIds,pfIds = gtsPFacet( 'Mesh_cylinder.gts',shift=(0,0,shiftfactor),scale=1,radius=rParticle,wire=False,fixed=False, materialNodes='gridNodeMat',material='pFacetMat',color=Vector3(0.5,1,0.5) ) # E.c.3). define different sections of membrane membrane_grid = [O.bodies[s] for s in nodesIds ] membrane_pfacet = [O.bodies[s] for s in pfIds] # E.c.4). run one interaction for f in membrane_grid: f.shape.color = Vector3(0,0,0) x = f.state.pos[0] y = f.state.pos[1] dist = math.sqrt(x*x+y*y) n = Vector3(x/dist,y/dist,0) z = f.state.pos[2] if z<=O.bodies[top_id].state.pos[2] and z>=O.bodies[bot_id].state.pos[2]: f.state.vel = -vel_ini_r*n O.engines[2].physDispatcher.functors[1].setCohesionNow=True while 1: O.run(1,True) break # E.c.5). redefine engine #**********************************************************************# O.engines=[ ForceResetter(), InsertionSortCollider([ Bo1_Sphere_Aabb(), Bo1_PFacet_Aabb(), Bo1_Facet_Aabb(), Bo1_GridConnection_Aabb() ]), InteractionLoop( [ Ig2_Sphere_Sphere_ScGeom6D(), Ig2_GridNode_GridNode_GridNodeGeom6D(), Ig2_GridConnection_GridConnection_GridCoGridCoGeom(), Ig2_Sphere_PFacet_ScGridCoGeom(), Ig2_Facet_Sphere_ScGeom6D() ], [ Ip2_FrictMat_FrictMat_FrictPhys(), Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label="cohesiveIp") ], [ Law2_ScGeom_FrictPhys_CundallStrack(), Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm =True,always_use_moment_law=False,label='cohesiveLaw'), Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(), Law2_ScGridCoGeom_FrictPhys_CundallStrack(), ],label="iloop" ), GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=25,timestepSafetyCoefficient=0.8), NewtonIntegrator(gravity=(0,0,0),damping=0.3,label='newton'), PyRunner(command='confining()',iterPeriod=1), PyRunner(command='plotAddData',iterPeriod=1) ] #**********************************************************************# # set final friction angle, enable cohesion setContactFriction(radians(finalFricDegree)) O.engines[2].lawDispatcher.functors[1].always_use_moment_law=True O.engines[2].physDispatcher.functors[1].setCohesionNow=True O.engines[2].physDispatcher.functors[1].setCohesionOnNewContacts=False # E.c.6). confining # some initial parameters p_a = 0 p_r = 0 e_ini = 0 V_ini = 0 R_avg = 0 H_ini = 0 porosity = 0 # velocity vel_ini_a = 0.05*rate*height_0 vel_ini_r = 0.05*rate*height_0 # loops (fast-slow) for reaching target confining stress while 1: O.run(10,True) (p_r,p_a,e_ini,V_ini,R_avg,H_ini,porosity)=confining() p_r = abs(p_r) pressure = max(p_r,p_a) dif = p_r-p_a if (p_a > isoStress/1000): vel_ini_a = -abs(vel_ini_a) if (p_r > isoStress/1000): vel_ini_r = -abs(vel_ini_r) else: vel_ini_r = abs(vel_ini_r) #elif (p_a <= isoStress/1000): if (p_r > isoStress/1000): vel_ini_a = abs(vel_ini_a) vel_ini_r = -abs(vel_ini_r) else: if (pressure<0.9*isoStress/1000): if dif > 5: vel_ini_a = 1.05*abs(vel_ini_a) #elif dif < -5: vel_ini_r = 1.05*abs(vel_ini_r) if (pressure>=0.85*isoStress/1000 and pressure<=isoStress/1000): if dif > 1: if dif > 5: vel_ini_a = 1.5*abs(vel_ini_a) else: vel_ini_a = 1.01*abs(vel_ini_a) elif dif < -1: if dif < -5: vel_ini_r = 1.5*abs(vel_ini_r) else: vel_ini_r = 1.01*abs(vel_ini_r) mean = (p_r+p_a)/2 unb=unbalancedForce() print 'p= ',p_r,' q= ',p_a,' porosity= ',porosity,' unbalanced force: ',unb if abs(isoStress/1000-mean)/(isoStress/1000)<0.15 and abs(dif) <15: print 'initial strain: ',e_ini print 'initial volume: ',V_ini print 'Confining stage I is finished!' #break while 1: for f in membrane_grid: f.state.blockedDOFs = 'xyzXYZ' O.run(1,True) (p_r,p_a,e_ini,V_ini,R_avg,H_ini,porosity)=confining() p_r = abs(p_r) pressure = max(p_r,p_a) dif = p_r-p_a if (p_a > isoStress/1000): vel_ini_a = -abs(vel_ini_a) if (p_r > isoStress/1000): vel_ini_r = -abs(vel_ini_r) else: vel_ini_r = abs(vel_ini_r) #elif (p_a <= isoStress/1000): if (p_r > isoStress/1000): vel_ini_a = abs(vel_ini_a) vel_ini_r = -abs(vel_ini_r) else: if (pressure<0.9*isoStress/1000): if dif > 5: vel_ini_a = 1.05*abs(vel_ini_a) #elif dif < -5: vel_ini_r = 1.05*abs(vel_ini_r) if (pressure>=0.9*isoStress/1000 and pressure<=isoStress/1000): if dif > 1: if dif > 5: vel_ini_a = 1.1*abs(vel_ini_a) else: vel_ini_a = 1.01*abs(vel_ini_a) #elif dif < -1: if dif < -5: vel_ini_r = 1.1*abs(vel_ini_r) else: vel_ini_r = 1.01*abs(vel_ini_r) mean = (p_r+p_a)/2 unb=unbalancedForce() print 'p= ',p_r,' q= ',p_a,' porosity= ',porosity,' unbalanced force: ',unb if abs(isoStress/1000-mean)/(isoStress/1000)<0.05 and abs(dif) <10: print 'initial strain: ',e_ini print 'initial volume: ',V_ini print 'Confining stage II is finished!' #break print 'V_ini= ',V_ini # E.c.7). stablize #**********************************************************************# O.engines=[ ForceResetter(), InsertionSortCollider([ Bo1_Sphere_Aabb(), Bo1_PFacet_Aabb(), Bo1_Facet_Aabb(), Bo1_GridConnection_Aabb() ]), InteractionLoop( [ Ig2_Sphere_Sphere_ScGeom6D(), Ig2_GridNode_GridNode_GridNodeGeom6D(), Ig2_GridConnection_GridConnection_GridCoGridCoGeom(), Ig2_Sphere_PFacet_ScGridCoGeom(), Ig2_Facet_Sphere_ScGeom6D() ], [ Ip2_FrictMat_FrictMat_FrictPhys(), Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label="cohesiveIp") ], [ Law2_ScGeom_FrictPhys_CundallStrack(), Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm =True,always_use_moment_law=False,label='cohesiveLaw'), Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(), Law2_ScGridCoGeom_FrictPhys_CundallStrack(), ],label="iloop" ), GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=25,timestepSafetyCoefficient=0.8), NewtonIntegrator(gravity=(0,0,0),damping=0.3,label='newton'), PyRunner(command='stable()',iterPeriod=1), PyRunner(command='plotAddData',iterPeriod=1) ] #**********************************************************************# O.engines[2].lawDispatcher.functors[1].always_use_moment_law=True O.engines[2].physDispatcher.functors[1].setCohesionNow=True O.engines[2].physDispatcher.functors[1].setCohesionOnNewContacts=False vel_a = abs(vel_ini_a) while 1: O.run(100,True) unb=unbalancedForce() print 'unbalanced force: ',unb e_ini = (top[0].state.displ()[2] - bot[0].state.displ()[2])/height_0 f1 = sum(O.forces.f(b)[2] for b in topcap) f2 = sum(O.forces.f(b)[2] for b in botcap) f4 = 0 r_cum = 0 count = 0 for b in membrane_grid: x = b.state.pos[0] y = b.state.pos[1] dist = math.sqrt(x*x+y*y) n = Vector3(x/dist,y/dist,0) a = 2*3.1416*dist/(360/alpha)*6*rParticle#*(1+e) z = b.state.pos[2] if z<=O.bodies[top_id].state.pos[2] and z>=O.bodies[bot_id].state.pos[2]: f_local = O.forces.f(b.id) length = math.sqrt(f_local[0]*f_local[0]+f_local[1]*f_local[1]+f_local[2]*f_local[2]) cos_theta = (n[0]*f_local[0]+n[1]*f_local[1]+n[2]*f_local[2])/length p_normal = (length*cos_theta/a) f4 += (p_normal) r_cum += dist count += 1 r_avg = r_cum/count-rParticle fa = abs(.5*(f2-f1)) p_r = f4/count/1000 p_a = fa/(3.1416*r_avg*r_avg)/1000 print 'pr=', p_r, ' pa=',p_a if unb <= thresholdvalue: #break ####################################################################### # F. APPLYING DEVIATOR LOADING # ####################################################################### # F.a). redifine engines O.engines=[ ForceResetter(), InsertionSortCollider([ Bo1_Sphere_Aabb(aabbEnlargeFactor=enlargefactor), Bo1_PFacet_Aabb(), Bo1_Facet_Aabb(), Bo1_GridConnection_Aabb() ]), InteractionLoop( [ Ig2_Sphere_Sphere_ScGeom6D(interactionDetectionFactor=enlargefactor), Ig2_GridNode_GridNode_GridNodeGeom6D(), Ig2_GridConnection_GridConnection_GridCoGridCoGeom(), Ig2_Sphere_PFacet_ScGridCoGeom(), Ig2_Facet_Sphere_ScGeom6D() ], [ Ip2_FrictMat_FrictMat_FrictPhys(), Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label="cohesiveIp") ], [ Law2_ScGeom_FrictPhys_CundallStrack(), Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm =True,always_use_moment_law=False,label='cohesiveLaw'), Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(), Law2_ScGridCoGeom_FrictPhys_CundallStrack(), ],label="iloop" ), GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=25,timestepSafetyCoefficient=0.8), NewtonIntegrator(gravity=(0,0,0),damping=0.3,label='newton') ] O.engines[2].lawDispatcher.functors[1].always_use_moment_law=True O.engines[2].physDispatcher.functors[1].setCohesionNow=True O.engines[2].physDispatcher.functors[1].setCohesionOnNewContacts=False O.engines=O.engines+[ PyRunner(command='addforce()',iterPeriod=1,label='force'), PyRunner(command='plotAddData()',iterPeriod=1,label='recorder'), PyRunner(command='checkrecord()',realPeriod=10,label='checker') ] # F.b). define the velocity of membrane walls to maintain the volume constant condition vel_a = final_rate*abs(vel_ini_a) vel_r = vel_a*.5*width/height Vel_r = vel_r conStress = p_r*1000 ####################################################################### # G. UTILITIES # ####################################################################### # G.a). time step (recommanded by YADE) O.dt=0.3*PWaveTimeStep() t = O.dt # G.b). funtion for plot plot.plots={'e':('q','p'),'e':('u','ev')} plot.plot() O.saveTmp() O.timingEnabled=1 def consolidation(Tv): #see your soil mechanics handbook... U=1 for k in range(50): M=pi/2*(2*k+1) U=U-2/M**2*exp(-M**2*Tv) return U triax.goal2=-11000 from yade import plot def history(): plot.addData(e22=-triax.strain[1]-zeroe22,e22_theory=drye22+(1-dryFraction)*consolidation((O.time-zeroTime)*Cv/wetHeight**2)*1000./modulus,t=O.time,p=flow.getPorePressure((0.5,0.1,0.5)),s22=-triax.stress(3)[1]-10000) from yade import plot plot.plots={'t':('e22','e22_theory',None,'s22','p')} plot.plot() O.saveTmp() O.timingEnabled=1 from yade import timing print "starting oedometer simulation" O.run(200,1) timing.stats() #################################################################### # G. RUN # #################################################################### print '=========================' print "start triaxial simulation" print '=========================' O.run() Cheers Sam

When	Who	Action	New State
2021-07-26	Robert Caulk	Request for more information	Needs information
2021-07-26	Sam	Give more information	Open
2021-07-26	Robert Caulk	Answer	Answered
2021-07-27	Sam	Reopen	Open
2021-07-27	Robert Caulk	Request for more information	Needs information
2021-08-12	Launchpad Janitor	Expire	Expired

Yade

History of question #698135