yade-users team mailing list archive

[Janek Kozicki <janek_listy@xxxxx>]

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
Description: PNG image

#!/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()