示例 4: 组合超级椭球的堆积模拟

该示例用于组合超级椭球的堆积模拟

组合超级椭球有形如真实颗粒(砂土)的如下显著的特性 (例如， 平直度，非对称性以及棱角度)。 组合超级椭球的形状函数可以通过组装八个超级椭球的形状参数得到: $$\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$$

::: {.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}$$ :::

其中 $r_x^+,r_y^+,r_z^+$ 和 $r_x^-,r_y^-,r_z^-$ 分别表示颗粒在正负主轴 x, y, z 的长度; $\epsilon_1, \epsilon_2$ 控制颗粒表面的棱角度, 对于凸颗粒来说，它们取值范围在 $(0,2)$。 图 Fig. 3.18{reference-type=“ref” reference=“figmultipolysuper”} 形象的描述了 $\epsilon_1$ 和 $\epsilon_2$ 是如何影响颗粒形状的。 Python包 _superquadrics_utils 提供了两个用来生成组合超级椭球的函数 (参见 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)

组合超级椭球及其形状参数 $r_x^+=1.0, r_x^-=0.5, r_y^+=0.8, r_y^- = 0.9, r_z^+ = 0.4, r_z^- = 0.6$ (a) $\epsilon_1=0.4, \epsilon_2=1.5$, (b) $\epsilon_1=\epsilon_2=1.0$, (c) $\epsilon_1=\epsilon_2=1.5$.

这里我们给出了组合超级椭球堆积模拟的示例脚本:

  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] #颗粒生成的盒子大小
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()]   # 接触力本构模型
   ),
   newton
]
O.dt=5e-5
#生成颗粒
Addlayer(a,eps,ap,p_mat,boxsize,num_s,num_t)#Note: these particles may intersect with each other when we generate them.

#显示设置。用户可以通过界面窗口来直接更改对应设置。
sudodem.qt.Gl1_PolySuperellipsoid.wire=False
sudodem.qt.Gl1_PolySuperellipsoid.slices=10
sudodem.qt.Gl1_Wall.div=0 #hide the walls

图 3.19{reference-type=“ref” reference=“figpolysuperpacking”} 给出了组合超级椭球堆积模拟的结果示意图。

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