Modelling flow around cylinder: particles fall into a closed cylinder

Asked by Daria on 2020-09-14

I have got case in which flow around cylinder is simulated. This is YADE-OpenFOAM coupling case. I have got closed cylinder, but at some time particles enter it and simulation fails. How can I fix it?

scriptYade.py
#####
from __future__ import print_function
import sys
from yadeimport import *
from yade.wrapper import *
from yade.utils import *

initMPI() #Initialize the mpi environment, always required.
fluidCoupling = yade.FoamCoupling(); #Initialize the engine
fluidCoupling.getRank(); #part of Initialization.

#example of spheres in shear flow : two-way point force coupling
class simulation():

 def __init__(self):
  O.periodic = True
  O.cell.setBox(1,0.4,0.2)

  numspheres=1000
  young = 1e10#5e6
  young_steel = 2e11
  density = 2#2000
  density_steel = 7700
  poisson = 0.3#0.5
  normalCohesion = 1e7
  shearCohesion = 1e7
  frictionAngle = radians(15)

  mat1 = CohFrictMat(normalCohesion=normalCohesion, shearCohesion=shearCohesion, young=young, poisson=poisson, frictionAngle=frictionAngle, density=density, momentRotationLaw=True, label='spheremat')
  O.materials.append(mat1)
  mat2 = CohFrictMat(normalCohesion=normalCohesion, shearCohesion=shearCohesion, young=young_steel, poisson=poisson, frictionAngle=frictionAngle, density=density_steel, momentRotationLaw=True, label='wallmat')
  O.materials.append(mat2)

  epsilon = 1e-08
  minval = 0 + epsilon
  maxx = 1.0 - epsilon
  maxy = 0.4 - epsilon
  maxz = 0.2 - epsilon
  #wall coords, use facets for wall BC:
  v0 = Vector3(minval, minval, minval)
  v1 = Vector3(minval,minval,maxz)
  v2 = Vector3(maxx,minval,minval)
  v3 = Vector3(maxx,minval,maxz)

  v4 = Vector3(minval,maxy,minval)
  v5 = Vector3(minval,maxy,maxz)
  v6 = Vector3(maxx,maxy,minval)
  v7 = Vector3(maxx, maxy, maxz)

  lf0 = facet(vertices=[v0,v1,v2], material='wallmat')
  O.bodies.append(lf0)
  lf1 = facet(vertices=[v1,v2,v3], material='wallmat')
  O.bodies.append(lf1)

  uf0 = facet(vertices=[v4,v5,v6], material='wallmat')
  O.bodies.append(uf0)
  uf1 = facet(vertices=[v5,v6,v7], material='wallmat')
  O.bodies.append(uf1)

  ff0 = facet(vertices=[v0,v2,v6], material='wallmat')
  O.bodies.append(ff0)
  ff1 = facet(vertices=[v0,v6,v4], material='wallmat')
  O.bodies.append(ff1)

  bf0 = facet(vertices=[v1,v3,v7], material='wallmat')
  O.bodies.append(bf0)
  bf1 = facet(vertices=[v1,v7,v5], material='wallmat')
  O.bodies.append(bf1)

  oriBody = Quaternion(Vector3(1,0,0),(0))
  radius = 0.05
  O.bodies.append(geom.facetCylinder((0.5,0.2,0.1),radius=radius,height=0.2,orientation=oriBody,segmentsNumber=100,wallMask=4,material='wallmat'))

  #spheres
  mn, mx= Vector3(minval + epsilon, minval + epsilon, minval + epsilon), Vector3(maxx / 2 - radius - epsilon, maxy - epsilon, maxz - epsilon)
  sp = pack.SpherePack();
  sp.makeCloud(mn,mx,rMean=0.003,rRelFuzz=0.5, num=numspheres)
  O.bodies.append([sphere(center,rad,material='spheremat') for center,rad in sp])
  mn, mx= Vector3(maxx / 2 + radius + epsilon, minval + epsilon, minval + epsilon), Vector3(maxx - epsilon, maxy - epsilon, maxz - epsilon)
  sp = pack.SpherePack();
  sp.makeCloud(mn,mx,rMean=0.003,rRelFuzz=0.5, num=numspheres)
  O.bodies.append([sphere(center,rad,material='spheremat') for center,rad in sp])
  sphereIDs = [b.id for b in O.bodies if type(b.shape)==Sphere]

  #coupling engine settings
  fluidCoupling.setNumParticles(len(sphereIDs))
  fluidCoupling.setIdList(sphereIDs)
  fluidCoupling.isGaussianInterp=False; #use pimpleFoamYade for gaussianInterp

  # Integrator
  newton=NewtonIntegrator(damping=0.0, gravity = (0.0 ,0.0, 0.0))
  # add small damping in case of stability issues.. ~ 0.1 max, also note : If gravity is needed, set it in constant/g dir.

  O.timingEnabled==True

  O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()], allowBiggerThanPeriod=True),
   InteractionLoop(
    [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Facet_Sphere_ScGeom6D()],
    [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=True)],
    [Law2_ScGeom6D_CohFrictPhys_CohesionMoment()],
   GlobalStiffnessTimeStepper(timestepSafetyCoefficient=0.3, label = "ts"),
   fluidCoupling, #to be called after timestepper
   PyRunner(command='sim.printMessage()', iterPeriod= 1, label='outputMessage'),
   VTKRecorder(fileName='yadep/3d-vtk-',recorders=['spheres', 'facets', 'intr', 'force'],virtPeriod=0.001)
  ]

 def printMessage(self):
  print("********************************YADE-ITER = " + str(O.iter) +" **********************************")
  print("********************************YADE-TIME = " + str(O.time) +" **********************************")

 def printState(self, n):
  print('blaaaaa ' + str(n))

 def printDt(self):
  print(O.dt)

 def setDt(self, dt):
  O.dt = dt

 def irun(self,num):
  O.run(num,1)

if __name__=="__main__":
 sim = simulation()
 sim.irun(10000000)
 fluidCoupling.killMPI()

import builtins
builtins.sim=sim

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Daria
Solved:
2020-09-17
Last query:
2020-09-17
Last reply:
Daria (romanovadi) said : #1

Radius of cylinder in YADE must be little bigger then in OpenFOAM, in other way there is area which is out of OpenFOAM calculation area in which particles can enter in YADE calculations. This results an error. This area exists because in YADE cylinder is made of planes.

##################################
oriBody = Quaternion(Vector3(1,0,0),(0))
radius = 0.051
O.bodies.append(geom.facetCylinder((0.5,0.2,0.1),radius=radius,height=0.2,orientation=oriBody,segmentsNumber=100,wallMask=4,material='wallmat'))