Large size cylinderConnection in periodic scene

Asked by Xin Li

Hi everyone,

I am facing some issues when using cylinderConnection with large size (radius > width of periodic cell) under periodic conditions.
1.Axial force is very high or 'Nan' when the length of segment >= 1/2 cell size (sometimes, other small length has the same issue). Is there a way to delete the contact geom and phys between nodes and gridconnections as I don't need the inner interactions of cylinder in my case.
2.Losing interactions with particles or some particles cross the cylinder and then interactions are detected after a period of time, this will cause explosion. Some specific number and length of segments of the cylinder can work normally (for example, in my MWE, length of each segment of the cylinder is cell width/2.1 and repeat 20 times). I have not tried the same size as cell, because the error 'float division by zero' from gridPfacet.py line140-i.phys.kn=E*math.pi*(radius**2)/L (two nodes located in the boundaries of different cells and their distance is zero?).
My MWE as below:
# encoding: utf-8

from builtins import range
import math
import glob, os
import numpy as np
from yade import pack, plot
from yade import geom
from yade import export
from yade import polyhedra_utils
from yade.gridpfacet import *

# Soil sphere parameters
E_m=8.1654644182
E=pow(10,E_m) # micro Young's modulus
v=1.1131826506e-01 # micro Poisson's ratio
mu=2.6666064103e+01 # interparticle friction angle
beta=7.2094747073e-01 # rolling stiffness
eta=6.9562665265e-02 # plastic limit of rolling
rho = 2650 # soil density

spMat = O.materials.append(CohFrictMat(young=E,poisson=v,frictionAngle=radians(mu),density=rho,isCohesive=False,alphaKr=beta,alphaKtw=beta,momentRotationLaw=True,etaRoll=eta,etaTwist=eta))
wallMat = O.materials.append(CohFrictMat(young=E,poisson=v,frictionAngle=radians(0),density=rho,isCohesive=False,alphaKr=0,alphaKtw=0,momentRotationLaw=False,etaRoll=0,etaTwist=0))
extcylMat = O.materials.append(CohFrictMat(young=E,poisson=v,frictionAngle=radians(0),density=1190,isCohesive=False,alphaKr=0,alphaKtw=0,momentRotationLaw=True,etaRoll=0,etaTwist=0))
intcylMat = O.materials.append(CohFrictMat(young=0,poisson=0,frictionAngle=radians(0),density=1190,isCohesive=False,alphaKr=0,alphaKtw=0,momentRotationLaw=True,etaRoll=0,etaTwist=0))

O.periodic = True
x=0.024 # 0.02
y=0.1
z=0.1
O.cell.setBox(x,y*3,z*3)

D=0.075 # 37.5
r=0.00225 # 0.0015
d=r*2
sp=pack.SpherePack()
mn,mx=Vector3(0,y,z+D),Vector3(x,y*2,z*2)
sp.makeCloud(mn,mx,rMean=r,rRelFuzz=0,seed=1)
sphereIds1=sp.toSimulation(material=spMat)
print(len(sp))

for i in range(0,3):
    O.bodies.append(utils.box((x/6+x*i/3,y,z*1.5),(x/6,0,z/2),wire=True,fixed=True,color=[1,0,0],material=wallMat))
    O.bodies.append(utils.box((x/6+x*i/3,y*2,z*1.5),(x/6,0,z/2),wire=True,fixed=True,color=[1,0,0],material=wallMat))
    O.bodies.append(utils.box((x/6+x*i/3,y*1.5,z),(x/6,y/2,0),wire=True,fixed=True,color=[1,0,0],material=wallMat))

#define engines
O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_GridConnection_Aabb()],allowBiggerThanPeriod=True),
   InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D(),
         Ig2_Sphere_GridConnection_ScGridCoGeom(),
        ],
  [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True, setCohesionOnNewContacts=False)],
  [Law2_ScGridCoGeom_CohFrictPhys_CundallStrack(),
         Law2_ScGeom6D_CohFrictPhys_CohesionMoment(
     always_use_moment_law=True,
     useIncrementalForm=True
  )],
   ),
   GlobalStiffnessTimeStepper(timestepSafetyCoefficient=0.8),
   NewtonIntegrator(gravity=(0,0,-9.81),damping=0.2,label='newton'),
   PyRunner(initRun=False,iterPeriod=20000,command='check()'),
]

nodesIds=[]
cylIds=[]
vertices=[]
L=x/2.1
for i in range(0,21):
    vertices.append([i*L,y*1.5,z+0.5*D])
cylinderConnection(vertices,radius=0.5*D,fixed = True,dynamic=False,nodesIds=nodesIds, cylIds=cylIds, intMaterial=intcylMat, extMaterial=extcylMat)
print(nodesIds,cylIds)

def check():
    for i in O.interactions:
        if i.id1 in nodesIds and i.id2 in nodesIds:
            print(i.id1,i.id2,i.phys.normalForce.norm(),i.phys.shearForce.norm())
        if i.id1 in cylIds and i.id2 in cylIds:
            print(i.id1,i.id2,i.phys.normalForce.norm(),i.phys.shearForce.norm())
    print('unb=',unbalancedForce())
    if unbalancedForce()<0.01:
        print('stable')
        O.pause()

Take my MWE for example,settings of cylinderConnection remain the same and I tried several tests as below:
1.cell width=0.024m, cylinder radius=0.075m, minimum particle radius=0.00225mm-work normally
2.cell width=0.024m, cylinder radius=0.075m, minimum particle radius=0.00150mm-explosion
3.cell width=0.024m, cylinder radius=0.0375m, minimum particle radius=0.00225mm-explosion
4.cell width=0.024m, cylinder radius=0.0375m, minimum particle radius=0.00150mm-explosion
5.cell width=0.02m, cylinder radius=0.075m, minimum particle radius=0.00225mm-work normally
6.cell width=0.02m, cylinder radius=0.075m, minimum particle radius=0.00150mm-explosion
7.cell width=0.02m, cylinder radius=0.0375m, minimum particle radius=0.00225mm-work normally
8.cell width=0.02m, cylinder radius=0.0375m, minimum particle radius=0.00150mm-work normally
I am confused about how to set the cylinder especially when I need larger cylinder and smaller particle.
My yade version is Yade 20220218-6333~d331682~focal1 and system is Ubuntu 20.04.3 LTS
Dose anyone has the same problem when using large size cylinderConnection under periodic condition?
Thanks for any help.
Relative questions: https://answers.launchpad.net/yade/+question/181411; https://answers.launchpad.net/yade/+question/689446

Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Bruno Chareyre (bruno-chareyre) said :
#1

Hi,
I've seen working examples where the cylinder was discretized into multiple connexions, all smaller than period. It could workaround your issue.
Even this trick is not tested extensively, yo uwill have to try and see by yourself.
Regards
BRuno

Revision history for this message
Xin Li (portal3) said :
#2

Hi Bruno,
Thank you for your reply! Yes, the way to discretize the cylinder into specific number of connections works in my case. Now, I find if there are two gridConnections in a cell and each one has the size of half of width of periodic cell and repeat manually in 3 cells, most of my models can run without loosing interactions. I will try more tests.
Regards,
Xin

Revision history for this message
Bruno Chareyre (bruno-chareyre) said :
#3

Welcome.
If you define at least 3 (or maybe 4) you could turn off "allowBiggerThanPeriod".
But as long it detects contacts properly it should be the same.
B

Revision history for this message
Xin Li (portal3) said :
#4

Hi Bruno
I may still need the flag “AllowBiggerThanPeriod” even the distance between each node below cell width/3, 4 or 6 or higher number as the bounding box of each connection has a larger size > half of the cell size. Unless I cut the radius part of the Aabb of gridConnection on the axial direction under periodic scene from grid.cpp.

O = scene->cell->unshearPt(O);
O2 = scene->cell->unshearPt(O2);
O2 = O2 + scene->cell->hSize * GC->cellDist.cast<Real>(); for (int k = 0; k < 3; k++) {
if (k<1) {aabb->min[k] = min(O[k], O2[k]); aabb->max[k] = max(O[k], O2[k])} else {
aabb->min[k] = min(O[k], O2[k]) - GC->radius;
aabb->max[k] = max(O[k], O2[k]) + GC->radius;
}

Thanks
Xin

Revision history for this message
Bruno Chareyre (bruno-chareyre) said :
#5

Indeed, enlarged bboxes and cell distorsions may make 4 not enough, but you get the idea.
B

Can you help with this problem?

Provide an answer of your own, or ask Xin Li for more information if necessary.

To post a message you must log in.