Example 4: packing of poly-superellipsoids

Here’s where your user finds out if your project is for them.

Poly-superellipsoid has a capability of capturing major features (e.g., flatness, asymmetry, and angularity) of realistic particles such as sands in nature. A poly-superellipsoid is constructed by assembling eight pieces of superellipoids, governed by the following surface function at a local Cartesian coordinate system: $$\Big(\big|\frac{x}{r_{x}}\big|^{\frac{2}{\epsilon_{1}}} + \big|\frac{y}{r_{y}}\big|^{\frac{2}{\epsilon_{1}}}\Big)^{\frac{\epsilon_{1}}{\epsilon_{2}}} +\big |\frac{z}{r_{z}}\big|^{\frac{2}{\epsilon_{2}}} = 1$$ with

::: {.subequations} $$\begin{aligned} r_{x} = r_x^+ \quad{}\text{if}\quad{} x\geq 0 \quad{}\text{else}\quad{} r_x^- \label{appa} \ r_{y} = r_y^+ \quad{}\text{if}\quad{} y\geq 0 \quad{}\text{else}\quad{} r_y^- \label{appb} \ r_{z} = r_z^+ \quad{}\text{if}\quad{} z\geq 0 \quad{}\text{else}\quad{} r_z^- \label{appc}\end{aligned}$$ :::

where $r_x^+,r_y^+,r_z^+$ and $r_x^-,r_y^-,r_z^-$ are the principal elongation along positive and negative directions of x, y, z axes, respectively; $\epsilon_1, \epsilon_2$ control the squareness or blockiness of particle surface, and their possible values are in $(0,2)$ for convex shapes. Fig. 3.18{reference-type=“ref” reference=“figmultipolysuper”} intuitively shows how $\epsilon_1$ and $\epsilon_2$ affect particle surface. The module _superquadrics_utils provides two functions to generate a poly-superellipsoid (see Sec. 5.1.2{reference-type=“ref” reference=“modsuperquadrics”}):

  • NewPolySuperellipsoid(eps, rxyz, mat, rotate, isSphere, inertiaScale = 1.0)

  • NewPolySuperellipsoid_rot(eps, rxyz, mat, $q_w$, $q_x$, $q_y$, $q_z$, isSphere, inertiaScale = 1.0)

Poly-superellipsoids with $r_x^+=1.0, r_x^-=0.5, r_y^+=0.8, r_y^- = 0.9, r_z^+ = 0.4, r_z^- = 0.6$ and (a) $\epsilon_1=0.4, \epsilon_2=1.5$, (b) $\epsilon_1=\epsilon_2=1.0$, (c) $\epsilon_1=\epsilon_2=1.5$.

Here we give an exemplified script for packing of poly-superellipsoids as follows:

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118

#
from sudodem._superquadrics_utils import *
import math
import random
import numpy as np

a=1.1734e-2
isSphere=False 
ap=0.4
eps=0.5
num_s=200
num_t=8000
trials=0


wallboxsize=[10e-1,0,0] #size of the layer for the particles generation
boxsize=np.array([7e-1,7e-1,20e-1])

def down():
    for b in O.bodies:
        if(isinstance(b.shape,PolySuperellipsoid)):
            if(len(b.intrs())==0):
                b.state.vel=[0,0,-2]
                b.state.angVel=[0,0,0]
    
def gettop():
    zmax=0
    for b in O.bodies:
        if(isinstance(b.shape,PolySuperellipsoid)):
            if(b.state.pos[2]>=zmax):
                zmax=b.state.pos[2]
    return zmax


def Genparticles(a,eps,ap,mat,boxsize,num_s,num_t,bottom):
    global trials
    gap=max(a,a*ap)
    #print(gap)
    num=0
    coor=[]
    width=(boxsize[0]-2.*gap)/2.
    length=(boxsize[1]-2.*gap)/2.
    height=(boxsize[2]-2.*gap)
    #iteration=0
    
    while num<num_s and trials<num_t:
        isOK=True
        pos=[0]*3
        pos[0]=random.uniform(-width,width)
        pos[1]=random.uniform(gap-length,length-gap)
        pos[2]=random.uniform(bottom+gap,bottom+gap+height) 
          
        for i in range(0,len(coor)):
            distance=sum([((coor[i][j]-pos[j])**2.) for j in range(0,3)])
            if(distance<(4.*gap*gap)):
                
                isOK=False
                break
        
        if(isOK==True):        
            coor.append(pos)
            num+=1
            trials+=1
            #print(num)
    return coor 
               
    
def Addlayer(a,eps,ap,mat,boxsize,num_s,num_t):
   bottom=gettop()   
   coor=Genparticles(a,eps,ap,mat,boxsize,num_s,num_t,bottom)
   for b in coor:
      #xyz = np.array([1.0,0.5,0.8,1.5,1.2,0.4])*0.1
      xyz = np.array([1.0,0.5,0.8,1.5,0.2,0.4])*0.1
      bb=NewPolySuperellipsoid([1.0,1.0],xyz,mat,True,isSphere)#
      bb.state.pos=[b[0],b[1],b[2]]
      bb.state.vel=[0,0,-1]
      O.bodies.append(bb)   
   down()



p_mat = PolySuperellipsoidMat(label="mat1",Kn=1e5,Ks=7e4,frictionAngle=math.atan(0.3),density=2650,betan=0,betas=0) #define Material with default values
wall_mat = PolySuperellipsoidMat(label="wallmat",Kn=1e6,Ks=7e5,frictionAngle=0.0,betan=0,betas=0) #define Material with default values
wallmat_b = PolySuperellipsoidMat(label="wallmat",Kn=1e6,Ks=7e5,frictionAngle=math.atan(1),betan=0,betas=0) #define Material with default values

O.materials.append(p_mat)
O.materials.append(wall_mat)
O.materials.append(wallmat_b)


O.bodies.append(utils.wall(-0.5*wallboxsize[0],axis=0,sense=1, material = wall_mat))#left x
O.bodies.append(utils.wall(0.5*wallboxsize[0],axis=0,sense=-1, material = wall_mat))#right x
O.bodies.append(utils.wall(-0.5*boxsize[1],axis=1,sense=1, material = wall_mat))#front y
O.bodies.append(utils.wall(0.5*boxsize[1],axis=1,sense=-1, material = wall_mat))#back y
	
O.bodies.append(utils.wall(0,axis=2,sense=1, material =wallmat_b))#bottom z
#O.bodies.append(utils.wall(boxsize[0],axis=2,sense=-1, material = wallmat_b))#top z

newton=NewtonIntegrator(damping = 0.3,gravity=(0.,0.,-9.8),label="newton",isSuperquadrics=2)#isSuperquadrics=2 for Poly-superellipsoids

O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_PolySuperellipsoid_Aabb(),Bo1_Wall_Aabb()],verletDist=0.2*0.1),
   InteractionLoop(
      [Ig2_Wall_PolySuperellipsoid_PolySuperellipsoidGeom(), Ig2_PolySuperellipsoid_PolySuperellipsoid_PolySuperellipsoidGeom()], 
      [Ip2_PolySuperellipsoidMat_PolySuperellipsoidMat_PolySuperellipsoidPhys()], # collision "physics"
      [PolySuperellipsoidLaw()]   # contact law -- apply forces
   ),
   newton
]
O.dt=5e-5
#adding particles
Addlayer(a,eps,ap,p_mat,boxsize,num_s,num_t)#Note: these particles may intersect with each other when we generate them.

#setting the Display properties. You can go to the Display panel to set them directly.
sudodem.qt.Gl1_PolySuperellipsoid.wire=False
sudodem.qt.Gl1_PolySuperellipsoid.slices=10
sudodem.qt.Gl1_Wall.div=0 #hide the walls

After running the script above, the user is expected to see a packing as shown in Fig. 3.19{reference-type=“ref” reference=“figpolysuperpacking”}.

Packing of poly-superellipsoids. {#figpolysuperpacking width=“10cm”}