periodic boundary in axial direction for particles in cylinder pipe

Asked by Son Tung Dang on 2020-06-26

Hi.
I would like to set up the periodic boundary condition for cylinder pipe. I read some test cases given in examples folder. But most of them are in a box. I try to set up a test case as,

#!/usr/bin/python
# -*- coding: utf-8 -*-
from yade import pack,ymport,export,geom,bodiesHandling,qt
import array as arr
import numpy as np
import math
import random

import os
import errno

dp = 1e3
O.periodic=True

# Add material
O.materials.append(FrictMat(young=10e9,poisson=.25,frictionAngle=0.5,density=dp,label='Par'))

# Parameters, which will be passed into facets creators
kwMeshes={'color':[1,1,0],'wire':True,'dynamic':False,'material':0}
oriBody = Quaternion(Vector3(1,0,0),(pi/2.0))
# Cylinder
O.bodies.append(geom.facetCylinder((0.0,0.0,0.0),radius=0.05,height=0.8,orientation=oriBody,segmentsNumber=10,wallMask=4,**kwMeshes))

# Pack
sp=pack.SpherePack()
sp.makeCloud(minCorner=(-0.05,-0.25,-0.05),maxCorner=(0.05,0.25,0.05),rMean=.006,rRelFuzz=.5,periodic=True)
predicate = pack.inCylinder((0,-0.2,0),(0,0.2,0),0.0486)
sp = pack.filterSpherePack(predicate,sp,returnSpherePack=True)

for c,r in sp:
 O.bodies.append(utils.sphere(c,r))
# Make a wall
O.bodies.append(wall((0,0,-0.5),axis=2))
# Compute the number of particles
n=-1
for b in O.bodies:
 if isinstance(b.shape,Sphere):
  n=n+1
print(n)
# Set the velocity for particles
for b in O.bodies:
 if isinstance(b.shape,Sphere):
  v1=random.uniform(-1.0,1.0)*0.1
  v2=random.uniform(0,1.0)*0.0001
  v3=random.uniform(-1.0,1.0)*0.1
  b.state.vel=(v1,v2,v3)
O.engines=[
 # SubdomainBalancer(colorize=True,initRun=True,iterPeriod=100),
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()],label='collider',allowBiggerThanPeriod=True),
 # Represent the geometry of a contact point between twoSpheres
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()],
 ),
 NewtonIntegrator(damping=.1,exactAsphericalRot=True,gravity=(0.0,0.0,0.0)),
 qt.SnapshotEngine(fileBase='3d-',iterPeriod=200,label='snapshot'),
 PyRunner(command='finish()',iterPeriod=200000)
]
# O.dt=5.e-6
O.dt=PWaveTimeStep()
O.run(1,True)
# we must open the view explicitly (limitation of the qt.SnapshotEngine)
qt.View()

# this function is called when the simulation is finished
def finish():
 # snapshot is label of qt.SnapshotEngine
 # the 'snapshots' attribute contains list of all saved files
 makeVideo(snapshot.snapshots,'3d.mpeg',fps=10,bps=10000)
 O.pause()
#O.run(10000,True)
#from yade import timing
#timing.stats()
#quit()

But it gave a strange result.
Could you please figure out what is wrong?
FYI, I used the yade 2018.02b, python 2.7.17 and ubuntu 18.04.
Thank you.

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Son Tung Dang
Solved:
2020-06-29
Last query:
2020-06-29
Last reply:
2020-06-29
Jan Stránský (honzik) said : #1

Hello,

why exactly you find the results strange (cylinder split in 4 pieces, ...)?
What exactly are actually "results" (GUI view, some saved data, ...)
What is the meaning of the wall?

I have tried the code, but after O.step(), the step takes too much time, do you have same observations?

cheers
Jan

Son Tung Dang (sontd1090) said : #2

Hi Jan.
Thank you for your comment.
I have a cylinder splitting in 4 pieces. And it is a strange result I observed.
All best,
Tung

Jan Stránský (honzik) said : #3

The split in 4 pieces is perfectly OK because for the visualization purposes, the true position is wrapped into the periodic cell. If you have cylinder around point (0,0,0), its negative coordinates are displayed in corresponding "image" corners.
It has no effect on the simulation, which (being periodic) still simulates the cylinder as a compact cylinder.

If you don't like it, shift the simulation by half size of the periodic cell (to have the cylinder "in the middle") or switch to another visualization tool.

Just to make it clear, splitting in 4 pieces or shifting the simulation has no effect on the simulation, it is just visualization and its adjustment.

cheers
Jan

Son Tung Dang (sontd1090) said : #4

I got it. Thanks!
Could you please tell me how can I modify the size of periodic cell?
In addition, It also have the same observation as you are after O.step(). Is that normal?
All best,
Tung

Jan Stránský (honzik) said : #5

> how can I modify the size of periodic cell?

O.cell.setBox(x,y,z) # [1]

but the cell size has no effect on this "4 pieces problem" if the center of cylinder is at (0,0,0).
My suggestion was to let the cell as it is, just "shift" the simulation: increase all x and y coordinates by a number to shift the cylinder "in the middle".

> In addition, It also have the same observation as you are after O.step(). Is that normal?

no, it is not... I have no idea about the reason, especially as the first iteration is OK.. I have checked b.state.pos, vel, angVel and everything looks OK..
(sometimes I have this issue if my simulation is unstable and the particle fly away to "infinity", but it is not this case)

cheers
Jan

[1] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.Cell.setBox

Son Tung Dang (sontd1090) said : #6

I got it.
Thank you for your time and effort.
All best,
Tung