New question #695989 on Yade: https://answers.launchpad.net/yade/+question/695989
Hello, I am modeling compaction of two materials "modeled as spheres" with two particle size distribution. Material 1 sphere sizes is much smaller than material 2 as seen in the code below. When I run this code for ~20 seconds, spheres pass out the facets. If I comment this line "spFilter2.toSimulation(color=(0.152, 0.368, 0.988),material=Mat2)", the code works perfectly!. I have asked similar question [1] and was told to reduce O.dt, and use large stiffness for the facets. I have tried using O.dt=.01*utils.PWaveTimeStep() and a young = 1e20 for Mat 3 (facet materials) and I ran the code for few minuets. The smaller size spheres still escape the facet. Can anyone please explain why this problem occurs and how to solve it? Thanks, Othman [1] https://answers.launchpad.net/yade/+question/695929 Note: I am using yadedaily to run the code below. ----------------------- from yade import pack, export import numpy as np Mat1=O.materials.append(FrictMat(young = 1e6, poisson = 0.45,frictionAngle = 0.00349066, density=2340,label='portlandite')) Mat2=O.materials.append(FrictMat(young = 7.2e10, poisson = 0.17,frictionAngle = 0.558, density=2650,label='sand')) Mat3=O.materials.append(FrictMat(young = 1e14, poisson = 0.1)) SG1=2.34 ##portlandite SG2=2.65 ##ASTM quartz sand ##cylinder dimensions radiuscyl=(500e-6/2) heightcyl=610e-6 ##center of cylinder cx=0 cy=0 cz=0 ##Initial cube dimensions ### mnx=cx-(radiuscyl*1.1) mny=cy-(radiuscyl*1.1) mnz=0 mxx=cx+(radiuscyl*1.1) mxy=cy+(radiuscyl*1.1) mxz=heightcyl*1.1 volume=(pi*radiuscyl**2)*heightcyl ############################ spheres ############################# sp1=pack.SpherePack() sizes1=1e-6*np.array([15,15.76,22.94,25]) #modified to avoid spheres escaping facet due to large size difference, dt issues passing1=[0,0.75,0.9,1] sp1.makeCloud((mnx,mny,mnz),(mxx,mxy,mxz),psdSizes=sizes1,psdCumm=passing1,num=5811) #### cylinder extraction pred=pack.inCylinder((cx,cy,cz),(cx,cy,heightcyl),radiuscyl) spFilter1=filterSpherePack(pred,sp1,Material=Mat1, returnSpherePack=True) print (len (spFilter1)) spFilter1.toSimulation(color=(0.533, 0.803, 0.780),material=Mat1) mass1=utils.getSpheresMass() #### sizes and distribution are from gradation curve of sand ####### sp2=pack.SpherePack() sizes2=1e-6*np.array([149,150,300,425,600,850]) #Diameters of portlandite passing2=[0,0.02,0.23,0.675,0.98,1] sp2.makeCloud((mnx,mny,mnz),(mxx,mxy,mxz),psdSizes=sizes2,psdCumm=passing2,num=51) #### cylinder extraction pred=pack.inCylinder((cx,cy,cz),(cx,cy,heightcyl),radiuscyl) spFilter2=filterSpherePack(pred,sp2,Material=Mat2, returnSpherePack=True) print (len (spFilter2)) spFilter2.toSimulation(color=(0.152, 0.368, 0.988),material=Mat2) mass2=utils.getSpheresMass()-mass1 print ('mass 1 mass2 in g',mass1*1000, mass2*1000) print ('mass2/mass1 ', mass2/mass1) total_mass=utils.getSpheresMass() ############################ facets ############################# facets=geom.facetCylinder((cx,cy,heightcyl/2),radiuscyl,heightcyl,segmentsNumber=150,wallMask=4,color=(0.156, 0.164, 0.152),material=Mat3) cylinder=O.bodies.append(facets) yade.qt.View() ##creating disks d1=geom.facetCylinder((cx,cy,heightcyl),radiuscyl*0.99,0,segmentsNumber=150,wallMask=1,color=(0.156, 0.164, 0.152),material=Mat3) d2=geom.facetCylinder((cy,cx,cz),radiuscyl*0.99,0,segmentsNumber=150,wallMask=1,color=(0.156, 0.164, 0.152),material=Mat3) disk1IDs= O.bodies.append(d1) disk2IDs= O.bodies.append(d2) for i in disk1IDs: body= O.bodies[i] body.state.vel = (0,0,-1) for n in disk2IDs: body= O.bodies[n] body.state.vel = (0,0,0) ############################# compaction ############################# O.dt=.1*utils.PWaveTimeStep() O.engines=[ ForceResetter(), InsertionSortCollider([ Bo1_Sphere_Aabb(), Bo1_Facet_Aabb() ]), InteractionLoop( [ Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom(), ], [ Ip2_FrictMat_FrictMat_FrictPhys(), Ip2_FrictMat_FrictMat_FrictPhys(), ], [ Law2_ScGeom_FrictPhys_CundallStrack(), ], ), NewtonIntegrator(damping=.3), PyRunner(iterPeriod=2500,command='force()',initRun=True), ] O.run() stress=[] def force(): f1= [O.forces.f(i)[2] for i in disk1IDs] f=np.mean(f1) s=f/(pi*radiuscyl**2) stress.append(s) thickness=(O.bodies[disk1IDs[1]].state.pos[2])-(O.bodies[disk2IDs[1]].state.pos[2]) packing_density=utils.getSpheresVolume()/(thickness*pi*radiuscyl**2) print("stress, thickness, packing density ",s,thickness,packing_density) -- 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
