← Back to team overview

yade-users team mailing list archive

[Question #429604]: Cylinder and periodic boundary conditions

 

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

Hi all, 

I am having a problem when I use a cylinder with periodic boundary condition. 

I made a simple script to reproduce the problem. The script is a column of grain falling inside a box under gravity, with a cylinder fixed on the lower (infinite) plane of the box. When I do that, everything works fine (option periodic = 0 and cylinderBottom = 1 in the script below). 
If now, I remove the sides of the box and replace it  with periodic boundary condition (periodic = 1 and cylinderBottom = 1 in script below), some particles literally goes inside the cylinder. And this whatever the stiffness of the interaction, the time step used, the contact law used (Law2_ScGeom_ViscElPhys_Basic or Law2_ScGeom_FrictPhys_CundallStrack) and the length of the cylinder (which I made much larger than the periodic cell).  
This does not happen if I consider instead of the cylinder a row of fixed spheres (periodic = 1 and cylinderBottom = 0 in the script below). 

There is also a strange behavior when defining the cylinder position and extents, which might be linked to the present problem:
- It is not possible to define a cylinder which extents corresponds to an integer of the periodic cell extent, e.g. cylinder(begin = (L/2.,-100*L,ground+dp/2.),end =(L/2.,100*L,ground+dp/2.),.... It considers that the length is zero at this moment, probably applying
However, if you interchange begin and end, it works, i.e. cylinder(begin = (L/2.,100*L,ground+dp/2.),end =(L/2.,-100*L,ground+dp/2.),....
- The 3D view shows very strange position of the cylinder, which do not correspond to its actual position, when you consider the position of the cylinder nodes.

Do you have any idea what the problem could be linked to  ? 

Raphael


Script:
from yade import pack
from yade.gridpfacet import *

dp= 6e-3

periodic = 1
cylinderBottom = 1

ground = 0.2*dp
H = 300*dp
L = 20*dp

O.engines = [
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_GridConnection_Aabb()],allowBiggerThanPeriod = True),
 InteractionLoop(
    [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom(),Ig2_Sphere_GridConnection_ScGridCoGeom()],
    [Ip2_ViscElMat_ViscElMat_ViscElPhys()],
    [Law2_ScGeom_ViscElPhys_Basic()]
#    [Ip2_FrictMat_FrictMat_FrictPhys()],
#    [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 PyRunner(command = 'check()', virtPeriod = 0.1),
 NewtonIntegrator(damping=0., gravity = (0,0,-9.81))
 ]

#Material creation
O.materials.append(ViscElMat(en=0.5, et=1., young=5e6, poisson=0.5, density=2500., frictionAngle=0.4, label='Mat'))
#O.materials.append(FrictMat(young = 5e6, poisson = 0.5, density=2500,frictionAngle=0.4, label='Mat'))  

#Periodic Cell or containing box
if periodic:
	O.periodic = True 
	O.cell.setBox(L,L,H)
	O.bodies.append(box(center= (L/2.,L/2.,ground),extents=(100,100,0),fixed=True,color = (0.,1.,0.),wire = True,material = 'Mat'))
else:
	O.bodies.append(box(center= (L/2.,L/2.,ground),extents=(100,100,0),wire = True,fixed=True,color = (0.,1.,0.),material = 'Mat'))#Made invisible to see below
	O.bodies.append(box(center= (0,L/2.,H/2.),extents=(0,L/2.,H/2.),wire = True,fixed=True,color = (0.,1.,0.),material = 'Mat'))#Made invisible to see
	O.bodies.append(box(center= (L,L/2.,H/2.),extents=(0,L/2.,H/2.),fixed=True,color = (0.,1.,0.),material = 'Mat'))
	O.bodies.append(box(center= (L/2.,0,H/2.),extents=(L/2.,0,H/2.),wire = True,fixed=True,color = (0.,1.,0.),material = 'Mat'))#Made invisible to see
	O.bodies.append(box(center= (L/2.,L,H/2.),extents=(L/2.,0,H/2.),fixed=True,color = (0.,1.,0.),material = 'Mat'))

#Bottom fixed cylinder or fixed particles
if cylinderBottom:
	n = len(O.bodies)
	cylinder(begin = (L/2.,100*L,ground+dp/2.),end =(L/2.,-100*L,ground+dp/2.),radius = dp/2.,fixed = True,color = (0,0,1),intMaterial = 'Mat',extMaterial = 'Mat')#Made invisible to see inside
else:
	for n in range(0,int(L/dp)+1):
		O.bodies.append(sphere(center = (L/2.,n*dp,ground+dp/2.),radius = dp/2.,fixed = True,color = (0,0,1),wire = True,material = 'Mat'))#Made invisible to see inside

#Particle cloud for gravity deposition
partCloud = pack.SpherePack()
partCloud.makeCloud(minCorner=(4*L/10.,4*L/10.,ground+dp),maxCorner=(6*L/10.,6*L/10.,H),rRelFuzz=0., rMean=dp/2., num = 2000)
partCloud.toSimulation(material='Mat')

O.saveTmp()
O.dt = 1e-5

#Function to check if the center of a particle is contained inside the cylinder or pseudo cylinder
def check():
	for b in O.bodies:
		if b.dynamic and b.state.pos[2]<ground+dp:
			posCylXZ = Vector3(L/2.,0.,ground+dp/2.)
			pRelX =  b.state.pos[0]%L - posCylXZ[0]
			pRelZ =  b.state.pos[2] - posCylXZ[2]
			if sqrt(pow(pRelX,2) + pow(pRelZ,2))<dp/2.:#If the particle center is closer than the radius (i.e. if it is completely inside!)
				delta  = dp/2.-sqrt(pow(pRelX,2) + pow(pRelZ,2))
				print('particle {0} is inside the cylinder from an amount of {1} diameter'.format(b.id,delta/dp))
				O.pause()


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