yade-users team mailing list archive

["Alexander Eulitz \[Eugen\]" <question229901@xxxxxxxxxxxxxxxxxxxxx>]

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

Hi,
we talked about this some time before, when I reported about the problem about "loosing spheres" in a simulation (I uploaded a video)
Here is an example (corresponding code below). A sphere is falling caused by gravitation on the bottom of facetBox.
It becomes apparent that contact forces (i.e. force of the sphere in z direction) are dependent on the position the sphere hits the box.
If the sphere hits the box on a place where there two facets adjoin contact force will be higher than on any other place.
What is the reason for this?


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

from yade import pack, plot, geom, utils, timing

import os
import sys
import time

utils.readParamsFromTable(nr='NonBatch', height=0.05, r_sp=0.011980, m_sp=0.00351, contactLaw='HertzMindlin', E_sp=190e9, P_sp=0.30, Rho_sp=0.41, E_fl=190e9, P_fl=0.3, Rho_fl=0.41, density_fl=7835, dmp=0.2)

from yade.params.table import *
fileid=nr+str(height)+contactLaw
if utils.runningInBatch(): fileid=O.tags['description']

##globale Variablen
Messfrequenz=200000 #in Hz
plotVirtPeriod=1.0/(Messfrequenz) #noch nicht eingebaut.

#Density of the material [kg/m³] 
volume_sp=r_sp**3*4/3*math.pi
density_sp=m_sp/volume_sp
print 'Kugel:',volume_sp,'m³', density_sp
print 'Boden:',density_fl,'kg/m³'
 
O.materials.append(FrictMat(young=E_fl, poisson=P_fl, frictionAngle=Rho_fl, density=density_fl, label='mat_ground'))
O.materials.append(FrictMat(young=E_sp, poisson=P_sp, frictionAngle=Rho_sp, density=density_sp, label='mat_spheres'))

sp_left=O.bodies.append(utils.sphere([-0.15,0,height+r_sp],r_sp,color=[1,1,0], material='mat_spheres'))
sp_middle=O.bodies.append(utils.sphere([0,0,height+r_sp],r_sp,color=[1,1,0], material='mat_spheres'))
sp_right=O.bodies.append(utils.sphere([0.15,0,height+r_sp],r_sp,color=[1,1,0], material='mat_spheres'))
#ground=O.bodies.append(utils.sphere([0,0,-r_sp*100],r_sp*100, fixed=True, color=[0,1,0], material='mat_ground'))
#pot=O.bodies.append(utils.geom.facetCylinder((0,0,height/2),0.054,height,orientation=Quaternion((0, 0, 0),0),segmentsNumber=20,wallMask=2, color=(0,1,0), fixed=True, material='mat_ground'))
pot=O.bodies.append(utils.geom.facetBox((0,0,height/2),(0.3,0.3,height/2),orientation=Quaternion((0, 0, 0),0),wallMask=16, color=(0,1,0), fixed=True, material='mat_ground'))

print 'Programmstart am ',time.strftime("%d.%m.%Y um %H:%M:%S Uhr") 
print "Batch-Parameter"
print O.tags['params']

# Engines entsprechend Kontaktgesetz auswählen (quick&dirty)
if contactLaw == 'HertzMindlin':
	O.engines=[
	   ForceResetter(),
	   InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
	   InteractionLoop(
		  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
		  [Ip2_FrictMat_FrictMat_MindlinPhys()],
		  [Law2_ScGeom_MindlinPhys_Mindlin()]

	   ),
	  NewtonIntegrator(damping=dmp, gravity=(0,0,-9.81)),
	 PyRunner(virtPeriod=plotVirtPeriod,command="myAddPlotData()", label='plotter'),
	 PyRunner(virtPeriod=0.05,command="printTime()")
	]
elif contactLaw == 'HertzLinearShear':
	O.engines=[
	   ForceResetter(),
	   InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
	   InteractionLoop(
		  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
		  [Ip2_FrictMat_FrictMat_MindlinPhys()],
		  [Law2_ScGeom_MindlinPhys_HertzWithLinearShear()]

	   ),
	  NewtonIntegrator(damping=dmp, gravity=(0,0,-9.81)),
	 PyRunner(virtPeriod=plotVirtPeriod,command="myAddPlotData()", label='plotter'),
	 PyRunner(virtPeriod=0.05,command="printTime()")
	]
elif contactLaw == 'Cundall':
	O.engines=[
	   ForceResetter(),
	   InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
	   InteractionLoop(
		  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
		  [Ip2_FrictMat_FrictMat_FrictPhys()],
		  [Law2_ScGeom_FrictPhys_CundallStrack()]

	   ),
	  NewtonIntegrator(damping=dmp, gravity=(0,0,-9.81)),
	 PyRunner(virtPeriod=plotVirtPeriod,command="myAddPlotData()", label='plotter'),
	 PyRunner(virtPeriod=0.05,command="printTime()")
	]

O.dt=1*utils.PWaveTimeStep()
print 'Timestep:',O.dt

def myAddPlotData():
	Fmess_l=utils.sumForces([0], (0,0,1))
	Fmess_m=utils.sumForces([1], (0,0,1))
	Fmess_r=utils.sumForces([2], (0,0,1))
	force=plot.addData(t=O.time, Fnormal_left=Fmess_l, Fnormal_middle=Fmess_m, Fnormal_right=Fmess_r)

def printTime():
	print O.time
	

plot.plots={'t':('Fnormal_left', 'Fnormal_middle', 'Fnormal_right')}


O.saveTmp()

plot.plot()
utils.waitIfBatch()



-- 
You received this question notification because you are a member of
yade-users, which is an answer contact for Yade.