Question #688400 on Yade changed: https://answers.launchpad.net/yade/+question/688400
Akm gave more information on the question: I will post the full script here. I will attach the drive links for the particles because I exported those particles using a gts surface. I did the sample preparation in a separate file using radial expansion. So the particle positions are here-kindly download those text files. ''https://drive.google.com/drive/folders/1M2rKs- vnR3uHC_ybv2vliYLhvO2g0RTm?usp=sharing'' The complete script is as follows: from yade import pack, qt, export, ymport, plot #Default parameters or from table readParamsFromTable(noTableOk=True, # unknownOk=True, intfactor = 1.2, #Max value of 1.2 is preferred as it can be reverted back to 1 after a step. Shear_Vel=0.1, S_length=0.060, S_area=0.060*0.030, Sig0=300e3, #300kN/m2 kn=600e3, #300kPa/mm plotperiod=500, topname='/home/akm/Documents/Yade/install/bin/Full-10k-Top.txt', botname='/home/akm/Documents/Yade/install/bin/Full-10k-Bot.txt', spname='/home/akm/Documents/Yade/install/bin/snapshots-trial/plot', vtkname='/home/akm/Documents/Yade/install/bin/vtk-trial/s', ) from yade.params.table import * F_init=Sig0*S_area #Initial Normal force #MATERIAL SPECIFICATION mat1=CpmMat( young = 32e9, density=2400, poisson = 0.30, frictionAngle = radians(35), epsCrackOnset = 5e-1, relDuctility = 1.1, sigmaT = 1.2e6, ) mat2=CpmMat( young = 5e9, density=1940, poisson = 0.30, frictionAngle = radians(35), epsCrackOnset = 5e-2, relDuctility = 1.0001, sigmaT = 1.2e1, ) O.materials.append(mat1) O.materials.append(mat2) #UPPER PARTICLE IMPORT: upbox_particles=ymport.text(topname,material=mat2, color=(1,0,0)) O.bodies.append(upbox_particles) #LOWER PARTICLE IMPORT: lowbox_particles=ymport.text(botname,material=mat1, color=(0,1,0)) O.bodies.append(lowbox_particles) ################### ##GEOMETRY CREATION : #################### O.materials.append(FrictMat(frictionAngle=0,density=0,label='walls')) #Dimension of Packing d1,d2=Vector3(.01,.01,.01),Vector3(.07,.0783,.04) # corners of the initial packing xmin=d1[0] xmax=d2[0] ymin=d1[1] ymax=d2[1] zmin=d1[2] zmax=d2[2] ht_lower=0.03 ht_upper=0.03 asp_ht=0.0083 #Shear plate for shearing in the +X direction P_Shear=utils.box((xmin,ymin+(ht_lower/2)-(ht_lower/20),zmin+(zmax-zmin)/2 ), (0,(ht_lower/2),(zmax-zmin)/2 ),fixed=True, wire=True, color=(0,0,1)) O.bodies.append(P_Shear) P_Shear.state.blockedDOFs = "xyzXYZ" #Surrounding setup geometry U_box_1=geom.facetBox((xmin+(xmax-xmin)/2, ymin+ht_lower+asp_ht+(asp_ht+ht_upper)/2,zmin+(zmax-zmin)/2 ), ((xmax-xmin)/2,(asp_ht+ht_upper)/2,(zmax-zmin)/2 ), wallMask=2) # +X Boundary constraints U_box_2=geom.facetBox((xmin+(xmax-xmin)/2, ymin+ht_lower+(asp_ht*2+ht_upper)/2,zmin+(zmax-zmin)/2 ), ((xmax-xmin)/2,(asp_ht*2+ht_upper)/2,(zmax-zmin)/2 ), wallMask=1) # -X Boundary constraints L_box=geom.facetBox((xmin+(xmax-xmin)/2, ymin+(asp_ht+ht_lower)/2,zmin+(zmax-zmin)/2 ), ((xmax-xmin)/2,(asp_ht+ht_lower)/2,(zmax-zmin)/2 ), wallMask=6) Cover_box=geom.facetBox(((xmin+xmax)/2,(ymin+ymax)/2,(zmin+zmax)/2 ),((xmin+xmax+0.06)/2,(ymin+ymax+0.02)/2,(zmin+zmax+0.02)/2 ) ,wallMask=63) O.bodies.append(U_box_1) O.bodies.append(U_box_2) O.bodies.append(L_box) O.bodies.append(Cover_box) idlist=[] idlist.append(P_Shear.id) for i in L_box: idlist.append(i.id) #ids for the translation motion of the box for i in O.bodies: if isinstance(i.shape,Sphere): if i.mat==mat1: i.state.blockedDOFs = "yzXYZ" #blocking the mat1 movement to ensure that there is no moment at the end point #Adding a plate at the top: global plate plate_pos=max([b.state.pos[1]+(b.shape.radius) for b in O.bodies if isinstance(b.shape,Sphere)]) plate=utils.box((xmin+(xmax-xmin)/2,plate_pos,zmin+(zmax-zmin)/2), ((xmax-xmin)/2 - xmin/60, 0, (zmax-zmin)/2),fixed=True, wire=False,color=(0,0,1)) O.bodies.append(plate) print('Plate added') plate.state.blockedDOFs = "xyzXYZ" O.engines=[ ForceResetter(), InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intfactor,label='bo1s'),Bo1_Facet_Aabb(),Bo1_Wall_Aabb(),Bo1_Box_Aabb()]), InteractionLoop( [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intfactor,label='ig2s'),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()], [Ip2_FrictMat_FrictMat_FrictPhys(), Ip2_CpmMat_CpmMat_CpmPhys(cohesiveThresholdIter=1)], [Law2_ScGeom_FrictPhys_CundallStrack(), Law2_ScGeom_CpmPhys_Cpm()] ), NewtonIntegrator(gravity=(0,-9.81,0),damping=0.5), GlobalStiffnessTimeStepper(active=True,timeStepUpdateInterval=1,timestepSafetyCoefficient=0.8), PyRunner(iterPeriod=1,command='addPlotData()'), ] O.step() bo1s.aabbEnlargeFactor=1 ig2s.interactionDetectionFactor=1 def addPlotData(): Sf = -O.forces.f(P_Shear.id)[0] #Negative sign since the force applied is in -X direction Hor_dspl = 1000*P_Shear.state.displ()[0] #Hor disp in mm Shear_stress = Sf/S_area Vert_dspl = 1000*plate.state.displ()[1] Force_Stif=kn*Vert_dspl*S_area #Force(N)=Stiffness(N/m2/mm)*Displacement(mm)*Area(m2) Force_Reqd=Force_Stif+F_init #Final force on plate = Force due to CNS + Intial force on the plate(SigmaO) Forceonplate=O.forces.f(plate.id)[1] if Forceonplate<Force_Reqd: plate.state.vel=(0,-2,0) elif Forceonplate>Force_Reqd: plate.state.vel=(0,1,0) Norm_stress_applied=(Forceonplate)/S_area Norm_stress_calculated=Force_Reqd/S_area yade.plot.addData( Hor_Disp = Hor_dspl, Shear_stress = Shear_stress, Ver_Disp = Vert_dspl, Norm_stress_applied = Norm_stress_applied, Norm_stress_calculated=Norm_stress_calculated ) # note the space in 'i ' so that it does not overwrite the 'i' entry plot.plots={'Hor_Disp':('Shear_stress'),'Hor_Disp ':('Norm_stress_applied','Norm_stress_calculated'),} plot.plot() #EXTRA ENGINES: O.engines=O.engines+[TranslationEngine(ids=idlist,translationAxis=[1,0,0], velocity=Shear_Vel)] O.engines=O.engines+[SnapshotEngine(fileBase=spname,iterPeriod=plotperiod)] #O.engines=O.engines+[VTKRecorder(fileName=vtkname,recorders=['cpm','all'],iterPeriod=500)] yade.qt.View() O.saveTmp() #####Thanks in advance guys##### -- 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
