yade-users team mailing list archive

[chiara modenese <c.modenese@xxxxxxxxx>]

>
> If you still want to compare ratcheting in Law2_ScGeom_FrictPhys_Basic, a
>> quick hack is to change "true"->"false" at line ElasticContactLaw.cpp:77
>> Please tell me the result. :)
>>
>
Hi Bruno,

I revised a little bit the script (I attached it again) I sent you
concerning RotationEngine (which applies the rolling) and now it works well.
I preferred to maintain the use of Rotation and StepDisplacer engines, I did
not have time to modify NewtonIntegrator and anyway both solutions should
work.

I tested it using Law2_ScGeom_FrictPhys_Basic. Avoiding granular ratcheting
(which is the default-fixed option in the code atm) no shear force is
present at the end of the loop (just run the script as it is). Instead, I
tried the case in which we want to consider granular ratcheting (changing
directly the code as you suggest above) and a net shear force remains at the
end (in the case I was running, I have got a value of 0.183), although the
ball comes back exactly at the same initial position. Meanwhile, I have also
printed out the values of the branches applied to the rotations in ScGeom,
and yes they differ, in the case with ratcheting, between the first and
second rotation. They are exactly the same avoiding ratcheting. Maybe I try
also with other contact laws and I let you know.

Chiara

#!/usr/local/bin/yade-trunk -x
# -*- coding: utf-8 -*-
# -*- encoding=utf-8 -*-

#__________________________________________________________________
# scrpt to demonstrate the effect of granular ratcheting
from yade import utils

#__________________________________________________________________
# geometry
r1,r2=0.007,0.007 # radii
p1,p2=[0,0,0],[r1+0.9*r2,0,0] # center positions 

#__________________________________________________________________
# material
young=600.0e4 
poisson=0.6 
density=2.60e3 
frictionAngle=89.0

# append geometry and material
O.materials.append(FrictMat(young=young,poisson=poisson,density=density,frictionAngle=frictionAngle))
O.bodies.append(utils.sphere(p1,r1,dynamic=False,wire=True)) # body id=0 
O.bodies.append(utils.sphere(p2,r2,dynamic=False,wire=True)) # body id=1

SFixed=O.bodies[0]
SMoving=O.bodies[1]

#__________________________________________________________________
# list of engines
O.engines=[
	StepDisplacer(subscribedBodies=[1],setVelocities=True,label='translation'),
	RotationEngine(subscribedBodies=[1],rotationAxis=[0,0,1],label='rotation'), # rotate around the centroid (if rotateAroundZero=False as it is by default)
	PeriodicPythonRunner(iterPeriod=1,command='displ()'),
	ForceResetter(),
	BoundDispatcher([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
	InsertionSortCollider(),
	InteractionGeometryDispatcher([Ig2_Sphere_Sphere_ScGeom()]), 
	InteractionPhysicsDispatcher([Ip2_FrictMat_FrictMat_FrictPhys()]), 
	LawDispatcher([Law2_ScGeom_FrictPhys_Basic(label='contact')]), 
	#NewtonIntegrator(),	
	PeriodicPythonRunner(iterPeriod=1,command='data()'),
]

#__________________________________________________________________
# time step
O.dt=.02*utils.PWaveTimeStep()
O.saveTmp('init')

#__________________________________________________________________
from yade import qt
qt.View()
qt.Controller()
from yade import plot

def displ():
	i=O.interactions[0,1]
	if O.iter==2:
		#O.bodies[1].state.vel=(-100,0,0)
		#O.bodies[1].state.angVel=(0,0,0)
		translation.deltaSe3=(-0.001*i.geom.normal,Quaternion.Identity) # apply displacement
		rotation.angularVelocity=0.0 
	elif O.iter==3: 
		#O.bodies[1].state.vel=(0,0,0)
		#O.bodies[1].state.angVel=(0,15,0)
		translation.deltaSe3=(Vector3(0,0,0),Quaternion.Identity) 
		rotation.angularVelocity=5000. # apply rotation
	elif O.iter==4:
		#O.bodies[1].state.vel=(100,0,0)
		#O.bodies[1].state.angVel=(0,0,0)
		translation.deltaSe3=(0.001*i.geom.normal,Quaternion.Identity) # apply displacement (come back to initial position)
		rotation.angularVelocity=0.0 
	elif O.iter==5:
		#O.bodies[1].state.vel=(0,0,0)
		#O.bodies[1].state.angVel=(0,-15,0)
		translation.deltaSe3=(Vector3(0,0,0),Quaternion.Identity) 
		rotation.angularVelocity=-5000. # apply rotation (come back to initial position)
	elif O.iter==6:
		#O.bodies[1].state.vel=(0,0,0)
		#O.bodies[1].state.angVel=(0,0,0)
		translation.deltaSe3=(Vector3(0,0,0),Quaternion.Identity)
		rotation.angularVelocity=0.0 

def data():
	i=O.interactions[0,1]
	if O.iter==2:		
		it=O.iter; un=i.geom.penetrationDepth; Fn=i.phys.normalForce.norm(); Fs=i.phys.shearForce.norm()
		print 'it = %i; un = %e; Fn = %e; Fs = %e' % (it, un, Fn, Fs), '\n\n--- plus dx ---'
	elif O.iter==3:
		it=O.iter; un=i.geom.penetrationDepth; Fn=i.phys.normalForce.norm(); Fs=i.phys.shearForce.norm()
		print 'it = %i; un = %e; Fn = %e; Fs = %e' % (it, un, Fn, Fs), '\n\n--- plus dteta ---'
	elif O.iter==4:
		it=O.iter; un=i.geom.penetrationDepth; Fn=i.phys.normalForce.norm(); Fs=i.phys.shearForce.norm()	
		print 'it = %i; un = %e; Fn = %e; Fs = %e' % (it, un, Fn, Fs), '\n\n--- minus dx ---'
	elif O.iter==5:
		it=O.iter; un=i.geom.penetrationDepth; Fn=i.phys.normalForce.norm(); Fs=i.phys.shearForce.norm()
		print 'it = %i; un = %e; Fn = %e; Fs = %e' % (it, un, Fn, Fs), '\n\n--- minus dteta ---'
	else: 
		it=O.iter; un=i.geom.penetrationDepth; Fn=i.phys.normalForce.norm(); Fs=i.phys.shearForce.norm()
		print 'it = %i;  un = %e; Fn = %e; Fs = %e' % (it, un, Fn, Fs), '\n'

O.run(7,True)