yade-dev team mailing list archive

[chiara modenese <c.modenese@xxxxxxxxx>]

Hi all,

 I need your help to understand if my observation (see what follows) is
wrong and eventually why it would be the case.

If you remember I was asking how the ids of the contact are assigned. The
answer was “random”, or better we can identify two main passages:


   1.

   collider
   2.

   geometry functor (it could be assymetric and so accept objects only in
   one specific order)

 The reason for which I was interested in all of this relies in the fact
that to me the way we now compute the incremental shear force appears
related to the ids of the particles.

 I took a very simple example (I attach the py script, I am using r2049 but
it should work also with the last one, I am just using the elastic contact
law). A particle is initially positioned in contact with a box. Then I apply
a normal force and a shear one in order to get a tangential reaction at
contact (that is why the normal force).

 FIRST CASE (you should get this just running the script as it is)

***ids bodies:

sphere = 0

box = 1

 ***ids contact:

id1 = 1 (box)

id2 = 0 (sphere)

 ***Shear force (type F8): negative (at least for a while, then the sign
obviously changes during the simulation)

 SECOND CASE (Here I swap the ids at contact. How? I “force” them in the
contact law, you can also do that just inverting the ids in it)

***ids bodies:

sphere = 0

box = 1

 ***ids contact:

id1 = 0 (sphere)

id2 = 1 (box)

 ***Shear force (type F8): positive! → if you go on with the iterations (the
script runs only few of them) you will see we are no more opposing the
motion since the shear force, this time, holds the wrong sign.

 You can repeat the test, is simple. The fact that the shear force changes
sign from one case to the other, is not correct. The shear force is
dependent on the relative shear velocity that in turn depends on the ids at
contact as we define them. To me this happens because we do not define a
local system.

Please correct me if you see I am wrong or please explain me why the shear
force changes sign (if it was really a global system as was already pointed
out, it should not happen).

 Thanks, Chiara

#!/usr/local/bin/yade-trunk -x
# -*- coding: utf-8 -*-
# -*- encoding=utf-8 -*-
##
## script to test the shear force in the elastic plastic contact law

O.initializers=[BoundDispatcher([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()])]

#####################
## *** 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()]
	),
	PeriodicPythonRunner(iterPeriod=1,command='runEngine()'),
	NewtonIntegrator(damping=0.0),
	PeriodicPythonRunner(iterPeriod=1,command='myAddPlotData()'),
]

######################
## *** MATERIAL *** ##
######################
mat=FrictMat(young=600.0e6,poisson=0.6,density=2.60e3,frictionAngle=25,label='Friction')
O.materials.append(mat)

######################
## *** GEOMETRY *** ##
######################
from yade import utils
sphere=utils.sphere([0.1,0,0],0.1,color=[0,1,0],dynamic=True,wire=True,material='Friction')
sphere.state.blockedDOFs=['rx','ry','rz'] #block particles rotations
box=utils.box(center=[-0.025,0,0],extents=[.025,1.0,1.0],color=[0,1,0],dynamic=False,wire=True,material='Friction')
O.bodies.append(sphere)# id=0
O.bodies.append(box)# id=1
print "\nsphere id = ", sphere.id
print "box id = ", box.id, "\n"

######################
## *** TIMESTEP *** ##
######################
O.dt=.2*utils.PWaveTimeStep()
O.saveTmp('Elastic')

from yade import qt
qt.View()
qt.Controller()


############################################
##### now the part pertaining to plots #####
############################################

from math import *
from yade import plot
plot.plots={'time':('fn',),'time_':('Fshear',),'time__':('un',),'time___':('vrel',)}

## this function is called by plotDataCollector
## it should add data with the labels that we will plot
## if a datum is not specified (but exists), it will be NaN and will not be plotted

###########################################
## *** FUNCTIONS for PERIODIC RUNNER *** ##
###########################################
def runEngine():
	ForceEngine(force=(-1e5,2e2,0),subscribedBodies=[0])()

def myAddPlotData():
	i=O.interactions[1,0]
	## store some numbers under some labels
	plot.addData(fn=i.phys.normalForce[0],us=sphere.state.pos[1],Fshear=i.phys.shearForce[1],vrel=box.state.vel[1]-sphere.state.vel[1],un=sphere.state.pos[0],time=O.time,time_=O.time,time__=O.time,time___=O.time)	

O.run(20,True);

## We will have:
## 1) data in graphs (if you call plot.plot())
## 2) data in file (if you call plot.saveGnuplot('/tmp/a')
## 3) data in memory as plot.data['step'], plot.data['fn'], plot.data['un'], etc. under the labels they were saved