There's definitely a bug hidden somewhere. And it might be not in plasticDissipation.
See attached picture: it is 140'000 iterations. Total energy is.... not going through the roof. But isn't conserved either. It's just conserved more or less. How I did it: 1. I used neverErase flag, because something goes wrong when interaction breaks. So do not delete it, but only reset Fn=0 and Fs=0. [Law2_ScGeom_FrictPhys_Basic(label='dry',traceEnergy=True,neverErase=True)] the neverErase flasg shouldn't change anything in the resulting graphs of total energy. But evidently it helps! So this is a BUG #1 here. 2. But that total energy is still wrong. It increases, albeit a little slower. And for some unknown reason using shearDisp helps. As if the mistake of using shearDisp is correcting another mistake present somewhere else: plasticDissipation += (-shearDisp - (1/currentContactPhysics->ks)*(trialForce-shearForce))//plastic disp. .dot(shearForce);//active force And this is BUG #2 in fact I suppose that BUG#1 == BUG#2 :) But where it is? -- Janek Kozicki http://janek.kozicki.pl/
<<attachment: scr_3225.png>>
#!/home/janek/bin/yy --threads=1 # -*- coding: utf-8 -*- # -*- encoding=utf-8 -*- # Script to monitor the energy of a system from yade import utils #----------------------------------------------------- # material #----------------------------------------------------- young=600.0e6 poisson=0.6 density=2.60e3 frictionAngle=radians(0) O.materials.append(FrictMat(young=young,poisson=poisson,density=density,frictionAngle=frictionAngle,label='mat')) #----------------------------------------------------- # geometry #----------------------------------------------------- # create a random packing in a box space from yade import pack sp=pack.SpherePack() mn=Vector3(0,0,0) mx=Vector3(1.0,1.0,1.0) ntot=10 R=0.1 std=0.0 sp.makeCloud(minCorner=mn,maxCorner=mx,rMean=R,rRelFuzz=std,num=ntot,periodic=False,porosity=0.85) for s in sp: O.bodies.append(utils.sphere(s[0],s[1])) # create some boundaries wires=True O.bodies.append(utils.box(center=[-0.1,0.5,0.5],extents=[.0025,1.0,1.0],dynamic=False,wire=wires,material='mat')) O.bodies.append(utils.box(center=[1.1,0.5,0.5],extents=[.0025,1.0,1.0],dynamic=False,wire=wires,material='mat')) O.bodies.append(utils.box(center=[1.1/2,-0.1,0.5],extents=[1.0,0.0025,1.0],dynamic=False,wire=wires,material='mat')) O.bodies.append(utils.box(center=[1.1/2,1.1,0.5],extents=[1.0,0.0025,1.0],dynamic=False,wire=wires,material='mat')) O.bodies.append(utils.box(center=[0.5,0.5,-0.1],extents=[1.0,1.0,0.0025],dynamic=False,wire=wires,material='mat')) O.bodies.append(utils.box(center=[0.5,0.5,1.1],extents=[1.0,1.0,0.0025],dynamic=False,wire=wires,material='mat')) #----------------------------------------------------- # initial condition (velocities) #----------------------------------------------------- v=7.0 for b in O.bodies: if b.id%2 == 0: b.state.vel=Vector3(-v,v,-v) # assign an initial velocity else: b.state.vel=Vector3(v,-v,v) # assign an initial velocity #----------------------------------------------------- # list of engines #----------------------------------------------------- O.engines=[ ForceResetter(), BoundDispatcher([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]), InsertionSortCollider(), InteractionDispatchers( [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()], [Ip2_FrictMat_FrictMat_FrictPhys()], [Law2_ScGeom_FrictPhys_Basic(label='dry',traceEnergy=True,neverErase=True)] ), NewtonIntegrator(damping=0.0), # *** NO DAMPING *** PeriodicPythonRunner(iterPeriod=10,command='monitoring()') ] #----------------------------------------------------- # time step #----------------------------------------------------- O.dt=.1*utils.PWaveTimeStep() O.saveTmp('init') from yade import qt qt.View() qt.Controller() #----------------------------------------------------- # plot some results #----------------------------------------------------- from math import * from yade import plot plot.plots={'t':('Ek','Eel','Slip','Etot'),'t_':('Slip')} def monitoring(): plot.addData(Ek=utils.kineticEnergy(),Eel=dry.elasticEnergy(),Slip=dry.plasticDissipation(),t=O.time,t_=O.time, Etot=utils.kineticEnergy()+dry.elasticEnergy()+dry.plasticDissipation()) O.run(10000,True); print "\n**** Now friction angle is set from 0 to 25 degrees. ****\n"; O.materials[0].frictionAngle=radians(25) O.run(130000,True); plot.plot()
_______________________________________________ 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