8.6. 扩展超椭球的堆积模拟

8.6.1. 扩展超椭球

扩展超椭球（Poly-superellipsoid）可以很好捕捉自然界中沙子等真实颗粒的主要特征(如长细比、非对称性和棱角度)的能力。 扩展超椭球由八块超椭球组合而成的, 在局部笛卡尔坐标系下形状控制方程如下：

\[\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\]
\[\begin{split}r_{x} &= r_x^+ \quad{}\text{if}\quad{} x\geq 0 \quad{}\text{else}\quad{} r_x^- \\ r_{y} &= r_y^+ \quad{}\text{if}\quad{} y\geq 0 \quad{}\text{else}\quad{} r_y^- \\ r_{z} &= r_z^+ \quad{}\text{if}\quad{} z\geq 0 \quad{}\text{else}\quad{} r_z^-\end{split}\]

其中 \(r_x^+,r_y^+,r_z^+\) 和 \(r_x^-,r_y^-,r_z^-\) 分别表示颗粒在正负主轴 x, y, z 的长度; \(\epsilon_1, \epsilon_2\) 控制颗粒表面的棱角度, 对于凸颗粒来说，它们取值范围在 $(0,2)$。图 8.18 形象的描述了 \(\epsilon_1\) 和 \(\epsilon_2\) 是如何影响颗粒形状的。 Python包 _superquadrics_utils 提供了两个用来生成扩展超椭球的函数:

  • 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)

../_images/multipolysuper1.png

图 8.18 扩展超椭球及其形状参数 \(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\).

8.6.2. 建模脚本

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

  1. 导入基础包
    1from sudodem._superquadrics_utils import *
2import math
3import random
4import numpy as np
  2. 模拟参数
     1a=1.1734e-2
 2isSphere=False
 3ap=0.4
 4eps=0.5
 5num_s=200
 6num_t=8000
 7trials=0
 8# 颗粒生成的盒子大小
 9wallboxsize=[10e-1,0,0]
10boxsize=np.array([7e-1,7e-1,20e-1])
  3. 定义辅助函数
     1# 给定所有颗粒向下的速度
 2def down():
 3    for b in O.bodies:
 4        if(isinstance(b.shape,PolySuperellipsoid)):
 5            if(len(b.intrs())==0):
 6                b.state.vel=[0,0,-2]
 7                b.state.angVel=[0,0,0]
 8
 9# 获得堆积体最高位置
10def gettop():
11    zmax=0
12    for b in O.bodies:
13        if(isinstance(b.shape,PolySuperellipsoid)):
14            if(b.state.pos[2]>=zmax):
15                zmax=b.state.pos[2]
16    return zmax
17
18# 生成颗粒
19def Genparticles(a,eps,ap,mat,boxsize,num_s,num_t,bottom):
20    global trials
21    gap=max(a,a*ap)
22    #print(gap)
23    num=0
24    coor=[]
25    width=(boxsize[0]-2.*gap)/2.
26    length=(boxsize[1]-2.*gap)/2.
27    height=(boxsize[2]-2.*gap)
28    #iteration=0
29    while num<num_s and trials<num_t:
30        isOK=True
31        pos=[0]*3
32        pos[0]=random.uniform(-width,width)
33        pos[1]=random.uniform(gap-length,length-gap)
34        pos[2]=random.uniform(bottom+gap,bottom+gap+height)
35
36        for i in range(0,len(coor)):
37            distance=sum([((coor[i][j]-pos[j])**2.) for j in range(0,3)])
38            if(distance<(4.*gap*gap)):
39
40                isOK=False
41                break
42
43        if(isOK==True):
44            coor.append(pos)
45            num+=1
46            trials+=1
47            #print(num)
48    return coor
49
50# 添加一层颗粒
51def Addlayer(a,eps,ap,mat,boxsize,num_s,num_t):
52bottom=gettop()
53coor=Genparticles(a,eps,ap,mat,boxsize,num_s,num_t,bottom)
54for b in coor:
55    #xyz = np.array([1.0,0.5,0.8,1.5,1.2,0.4])*0.1
56    xyz = np.array([1.0,0.5,0.8,1.5,0.2,0.4])*0.1
57    bb=NewPolySuperellipsoid([1.0,1.0],xyz,mat,True,isSphere)#
58    bb.state.pos=[b[0],b[1],b[2]]
59    bb.state.vel=[0,0,-1]
60    O.bodies.append(bb)
61down()
  4. 定义材料
    1p_mat = PolySuperellipsoidMat(label="mat1",Kn=1e5,Ks=7e4,frictionAngle=math.atan(0.3),density=2650,betan=0,betas=0) #define Material with default values
2wall_mat = PolySuperellipsoidMat(label="wallmat",Kn=1e6,Ks=7e5,frictionAngle=0.0,betan=0,betas=0) #define Material with default values
3wallmat_b = PolySuperellipsoidMat(label="wallmat",Kn=1e6,Ks=7e5,frictionAngle=math.atan(1),betan=0,betas=0) #define Material with default values
4
5O.materials.append(p_mat)
6O.materials.append(wall_mat)
7O.materials.append(wallmat_b)
  5. 生成墙/容器
    1O.bodies.append(utils.wall(-0.5*wallboxsize[0],axis=0,sense=1, material = wall_mat))#left x
2O.bodies.append(utils.wall(0.5*wallboxsize[0],axis=0,sense=-1, material = wall_mat))#right x
3O.bodies.append(utils.wall(-0.5*boxsize[1],axis=1,sense=1, material = wall_mat))#front y
4O.bodies.append(utils.wall(0.5*boxsize[1],axis=1,sense=-1, material = wall_mat))#back y
5
6O.bodies.append(utils.wall(0,axis=2,sense=1, material =wallmat_b))#bottom z
7#O.bodies.append(utils.wall(boxsize[0],axis=2,sense=-1, material = wallmat_b))#top z
  6. 设置引擎
     1newton=NewtonIntegrator(damping = 0.3,gravity=(0.,0.,-9.8),label="newton",isSuperquadrics=2)#isSuperquadrics=2 for Poly-superellipsoids
 2
 3O.engines=[
 4ForceResetter(),
 5InsertionSortCollider([Bo1_PolySuperellipsoid_Aabb(),Bo1_Wall_Aabb()],verletDist=0.2*0.1),
 6InteractionLoop(
 7    [Ig2_Wall_PolySuperellipsoid_PolySuperellipsoidGeom(), Ig2_PolySuperellipsoid_PolySuperellipsoid_PolySuperellipsoidGeom()],
 8    [Ip2_PolySuperellipsoidMat_PolySuperellipsoidMat_PolySuperellipsoidPhys()], # collision "physics"
 9    [PolySuperellipsoidLaw()]   # 接触力本构模型
10),
11newton
12]
13O.dt=5e-5
  7. 生成颗粒
    1#Note: these particles may intersect with each other when we generate them.
2Addlayer(a,eps,ap,p_mat,boxsize,num_s,num_t)
  8. 显示设置

    用户可以通过界面窗口来直接更改对应设置。

    1sudodem.qt.Gl1_PolySuperellipsoid.wire=False
2sudodem.qt.Gl1_PolySuperellipsoid.slices=10
3sudodem.qt.Gl1_Wall.div=0 #hide the walls

图 8.19 给出了扩展超椭球堆积模拟的结果示意图。

../_images/packingpolysuper.png

图 8.19 超椭球堆积状态