New question #706768 on Yade: https://answers.launchpad.net/yade/+question/706768
Hi, I'm doing uniaxial compression of particles in cylinder container. I'm comparing 2 contact models (Hertz Mindlin and Cundall Strack). I noticed that the load needed to obtain volume ratio (volume of particles/volume of container) 1 is much higher when using Cundall Strack model (12 000 N, comparing to 2000 N when using Hertz Mindlin). All the parameters are the same in both simulations. Does anyone know why this happens? Also the load in both cases increases linearly, and by experiments after some volume ratio it should increase kind of exponentionally (for example, the steep of the load curve between volume ratios 0.85 and 1 is way steeper than for less volume ratios). Thanks in advance MWE: from yade import pack, plot, export import math #sphere parameters r=0.0025 rr=0.05 n=800 #material definition E=31e6 v=0.3897 Ro=1140 dmp=0.05 fr=15 O.materials.append(FrictMat(young=E, poisson=v, frictionAngle=radians(fr), density=Ro, label='spheres')) Ec=2.1e9 vc=0.33 Roc=7850 frc=22 O.materials.append(FrictMat(young=Ec, poisson=vc, frictionAngle=radians(frc), density=Roc, label='walls')) #cylinder container R=0.015 h=10*R c=(0,0,h/2) cylinder = geom.facetCylinder(center=c, radius=R, height=h, segmentsNumber=12, wallMask=6, material='walls') O.bodies.append(cylinder) #funnel definition c1=(0,0,h) dB=8*R dO=2*R hB=1.5*h hO=0.5*h hP=0 funnel = geom.facetBunker(center=c1, dBunker=dB, dOutput=dO, hBunker=hB,hOutput=hO, hPipe=hP, segmentsNumber=20, wallMask=4, material='walls') O.bodies.append(funnel) #sphere packing sp = pack.SpherePack() sp.makeCloud((-dB/3,-dB/3,1.5*h),(dB/3,dB/3,1.9*h), rMean = r, rRelFuzz = rr, num = n) sp.toSimulation(material = 'spheres') #setting simulation engines O.engines = [ ForceResetter(), InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb(), Bo1_Wall_Aabb()]), InteractionLoop( [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom(), Ig2_Wall_Sphere_ScGeom()], [Ip2_FrictMat_FrictMat_MindlinPhys()], [Law2_ScGeom_MindlinPhys_Mindlin(label='Mindlin')] ), NewtonIntegrator(gravity=(0, 0, -9.81), damping=dmp), PyRunner(command='checkUnbalanced()', realPeriod=2, label='checker'), ] #timestep O.dt = 0.3*PWaveTimeStep() #checking the status of gravity deposition and loading if done def checkUnbalanced(): if O.iter < 150000: return if unbalancedForce()>1: return #adding the loading plate O.bodies.append(wall(max([b.state.pos[2] + b.shape.radius for b in O.bodies if isinstance(b.shape, Sphere)]), axis=2, sense=-1)) global plate plate = O.bodies[-1] w=plate.state.pos[2] plate.state.vel =(0,0,-0.05) O.engines = O.engines + [PyRunner(command='addPlotData()', iterPeriod=100)] checker.command = 'unloadPlate()' #checking the loading status and unloading if done def unloadPlate(): V1=((2*R)*(2*R)*3.14/4)*plate.state.pos[2] V2=n*4*3.14*r*r*r/3 Vr=V2/V1 Fz = O.forces.f(plate.id)[2] S=Fz/(4*R*R*3.14/4) if Vr>=1.1: O.pause() #plotting the data def addPlotData(): if not isinstance(O.bodies[-1].shape, Wall): plot.addData() return Fz = O.forces.f(plate.id)[2] w=plate.state.pos[2] V1=((2*R)*(2*R)*3.14/4)*w V2=n*4*3.14*r*r*r/3 Vr=V2/V1 S=Fz/(4*r*r*3.14/4) plot.addData(Fz=Fz, w=plate.state.pos[2], unbalanced=unbalancedForce(), Vr=V2/V1, S=Fz/(4*R*R*3.14/4)) plot.plots = {'w':('Fz'),'Vr':('S'),} plot.plot() O.run() -- 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