yade-users team mailing list archive

[Steve <question668664@xxxxxxxxxxxxxxxxxxxxx>]

New question #668664 on Yade:
https://answers.launchpad.net/yade/+question/668664

Hi,

I am trying to insert a big clump into a bed of spheres, displace this clump in the spheres and record force and torque on this clump. 
However, the torque I get does not seem to make sense. 

I tried both O.forces.t and torque recorder.
The O.forces.t gives torque that does not match the force on the clump (the force is about 100 N and 0.15 meter away from the center of the clump but the torque recorded is only about 5N*m.

The torque recorder gives values that does not make much sense. I guess the value it recorded depends on the  rotation axis and zeroPoint. How are these defined w.r.t to? Are they in the global coordinates or the clump's local coordinates? 

Below is the code if you want to take a look at it.

Thanks,
Yifan


#!/usr/bin/python
# -*- coding: utf-8 -*-

import random
from yade import ymport
from yade import qt
from yade import plot
import sys,time

factor=1
test='2'
level=4700*factor
## PhysicalParameters 
sand=O.materials.append(FrictMat(young=3e8,poisson=0.3,density=1631,frictionAngle=0.5236))
#default=O.materials.append(FrictMat(young=1e7,poisson=.3,density=1e3,frictionAngle=.5))
default=O.materials.append(CohFrictMat(young=1e7,poisson=.3,density=1e3,frictionAngle=.5,alphaKr=2.0,alphaKtw=1.0,etaRoll=-1,isCohesive=False,momentRotationLaw=True))
# Clump
clump1=O.bodies.appendClumped([\
sphere([0,0,0.1], radius=0.1,material=default),\
sphere([0,0,0.15],radius=0.1,material=default),\
sphere([0,0,0.2], radius=0.1,material=default),\
sphere([0,0,0.25], radius=0.1,material=default),\
sphere([0,0,0.3], radius=0.1,material=default),\
sphere([0,0,0.35], radius=0.1,material=default),\
sphere([0,0,0.4], radius=0.1,material=default),\
])


#s=O.bodies.append(sphere(center=(0,0,0.5),radius=0.2,fixed=False))
#get clump ids:
ClumpID = clump1[0]
SpheresID = clump1[1]
O.bodies[ClumpID].state.blockedDOFs='XYZ'
O.bodies[ClumpID].state.angVel=(0,0,0)
O.bodies[ClumpID].state.vel=(0,0,0)
#print id_clump2
#Default Material:
#FrictMat(density=1e3,young=1e7,poisson=.3,frictionAngle=.5).

# Spheres
sphereRadius = 0.01
#nbSpheres = (45,45,60)
nbSpheres=(32,32,30)
# with 32,32,11 spheres settle around -0.08555m after 2s 

surfacesphere=[]
print "Creating %d spheres..."%(nbSpheres[0]*nbSpheres[1]*nbSpheres[2]),
for i in xrange(nbSpheres[0]):
    for j in xrange(nbSpheres[1]):
		for k in xrange(nbSpheres[2]):
			x = (i*2 - nbSpheres[0])*sphereRadius*1.1+sphereRadius*random.uniform(-0.1,0.1)
			y = (j*2 - nbSpheres[1])*sphereRadius*1.1+sphereRadius*random.uniform(-0.1,0.1)
			z = -k*sphereRadius*2.2-0.01
			r = random.uniform(sphereRadius,sphereRadius*0.9)
			fixed = False
			color=[0.51,0.52,0.4]
			if (i==0 or i==nbSpheres[0]-1 or j==nbSpheres[1]-1 or j==0 or k==nbSpheres[2]-1):
				fixed = True
				color=[0.21,0.22,0.1]
			tmp=O.bodies.append(sphere([x,y,z],r,color=color,fixed=fixed,material=default))
			if k==0 and not fixed:
				surfacesphere.append(tmp)

			
print "done\n"


## Estimate time step
#O.dt=PWaveTimeStep()
O.dt=0.0001



	
def changeforce():
	if O.iter <15000*factor:
		O.bodies[ClumpID].state.pos=(-0.1,0,0.25)
		O.bodies[ClumpID].state.angVel=(0,0,0)
		O.bodies[ClumpID].state.vel=(0,0,0)
	elif O.iter == 15000*factor:
		O.bodies[ClumpID].state.pos=(-0.1,0,0.25-0.238)
		O.bodies[ClumpID].state.angVel=(0,0,0)
		O.bodies[ClumpID].state.vel=(0,0,0)
	elif O.iter > 15000*factor and O.iter < 15000*factor+16000:
		O.bodies[ClumpID].state.pos=(-0.1,0,0.25-0.238-O.dt*0.1*(O.iter-15000))
	elif O.iter >= 15000*factor+16000:
		O.bodies[ClumpID].state.pos=(-0.1+0.1414*O.dt*(O.iter-31000),0,0.25-0.238-0.16)

	
	

def addPlotdata():
	plot.addData(i=O.iter,x=O.bodies[ClumpID].state.pos[0],z=O.bodies[ClumpID].state.pos[2],f=O.forces.f(ClumpID,sync=True),t=O.forces.t(ClumpID,sync=True))
	#print O.bodies[ClumpID].state.pos[2]

def savedata():
	#pass
	if O.iter > 40000*factor:
		plot.saveDataTxt('displacement'+test+'.dat')

O.engines=[
	## Resets forces and momenta the act on bodies
	ForceResetter(),
	
	##display all kinds of info here
	
	## Using bounding boxes find possible body collisions.
	InsertionSortCollider([
		Bo1_Sphere_Aabb(),
		Bo1_Facet_Aabb(),
	]),
	InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
		#[Ip2_FrictMat_FrictMat_FrictPhys()],
		[Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
		[Law2_ScGeom_FrictPhys_CundallStrack()]
	),
	NewtonIntegrator(damping=0.3,gravity=[0,0,-9.81]),
	
	## Save force onClump
	ForceRecorder(ids=[ClumpID],file='force'+test+'.dat',iterPeriod=1),
	TorqueRecorder(ids=[ClumpID],file='itorque'+test+'.dat',iterPeriod=1,rotationAxis=(0,1,0),zeroPoint=(0,0,1)),
	TorqueRecorder(ids=[ClumpID],file='iitorque'+test+'.dat',iterPeriod=1,rotationAxis=(0,1,0),zeroPoint=(0,0,0)),
	TorqueRecorder(ids=[ClumpID],file='iiitorque'+test+'.dat',iterPeriod=1,rotationAxis=(0,1,0),zeroPoint=(0,0,-1)),
	PyRunner(command='changeforce()',iterPeriod = 1),
	PyRunner(command='addPlotdata()',iterPeriod = 1),
	PyRunner(command='savedata()',iterPeriod = 100)	
]

nbIter=50000*factor

qt.View()

O.stopAtIter=nbIter
O.run()




-- 
You received this question notification because your team yade-users is
an answer contact for Yade.