8.5. GJK颗粒的堆积模拟

8.5.1. GJK颗粒

在 SudoDEM 中, 有五种形状——球体、多面体、圆锥体、圆柱体、和立方体可采用接触检测GJK算法，因此我们暂且使用名字 GJKparticle 对这些颗粒形状进行归类分组。在Python模块 _gjkparticle_utils ,提供了GJK颗粒的构造函数, 包括:

• GJKSphere(radius,*margin*,*mat*)

• GJKPolyhedron(vertices,*extent*,*margin*,*mat*, rotate)

• GJKCone(radius, height, margin, mat, rotate)

• GJKCylinder(radius, height, margin, mat, rotate)

• GJKCuboid(extent, margin,*mat*, rotate)

注: 为了保证计算效率与计算稳定性的平衡性, 我们引入 margin参数来放大颗粒。该参数应该足够大, 以确保碰撞可以在被准确检测，同时该参数也应足够小，以保证颗粒形状不失真。

8.5.2. 脚本建模

与 运行一个简单的模拟 类似, 我们将会模拟GJK颗粒在立方体盒子中的自由落体堆积过程。

进行模拟前，我们需要先导入如下包:

from sudodem import plot, _gjkparticle_utils
from sudodem import *
from sudodem._gjkparticle_utils import *
import random as rand
import math

定义生成颗粒的参数信息：

cube_v0 = [[0,0,0],[1,0,0],[1,1,0],[0,1,0],[0,0,1],[1,0,1],[1,1,1],[0,1,1]]
cube_v = [[i*0.01 for i in j] for j in cube_v0]
isSphere=False
num_x = 5
num_y = 5
num_z = 20
R = 0.01

定义颗粒生成栅格：

#R: 相邻节点的距离
#num_x: x轴节点数量
#num_y: y轴节点数量
#num_z: z轴节点数量
#返回所有节点的位置信息
def GridInitial(R,num_x=10,num_y=10,num_z=20):
   pos = list()
   for i in range(num_x):
      for j in range(num_y):
         for k in range(num_z):
            x = i*R*2.0
            y = j*R*2.0
            z = k*R*2.0
            pos.append([x,y,z])
   return pos

设置材料属性

# 颗粒材料 p_mat
p_mat = GJKParticleMat(label="mat1",Kn=1e4,Ks=7e3,frictionAngle=math.atan(0.5),density=2650,betan=0,betas=0)
# 墙材料
wall_mat = GJKParticleMat(label="wallmat",Kn=1e4,Ks=7e3,frictionAngle=0.0,betan=0,betas=0)
# 墙材料b
wallmat_b = GJKParticleMat(label="wallmat",Kn=1e4,Ks=7e3,frictionAngle=math.atan(1),betan=0,betas=0)
O.materials.append(p_mat)
O.materials.append(wall_mat)
O.materials.append(wallmat_b)

生成立方体盒子

# create the box and add it to O.bodies
#left wall along x axis
O.bodies.append(utils.wall(-R,axis=0,sense=1, material = wall_mat))
#right wall along x axis
O.bodies.append(utils.wall(2.0*R*num_x-R,axis=0,sense=-1, material = wall_mat))
#front wall along y axis
O.bodies.append(utils.wall(-R,axis=1,sense=1, material = wall_mat))
#back wall along y axis
O.bodies.append(utils.wall(2.0*R*num_y-R,axis=1,sense=-1, material = wall_mat))
#bottom wall along z axis
O.bodies.append(utils.wall(-R,axis=2,sense=1, material =wallmat_b))

生成模拟引擎并设置固定时步：

newton=NewtonIntegrator(damping = 0.3,gravity=(0.,0.,-9.8),label="newton",isSuperquadrics=4)

O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_GJKParticle_Aabb(),Bo1_Wall_Aabb()],verletDist=0.2*0.01),
   InteractionLoop(
      [Ig2_Wall_GJKParticle_GJKParticleGeom(), Ig2_GJKParticle_GJKParticle_GJKParticleGeom()],
      [Ip2_GJKParticleMat_GJKParticleMat_GJKParticlePhys()], # collision "physics"
      [GJKParticleLaw()]   # contact law
   ),
   newton
]
O.dt=1e-5

多面体颗粒生成如下所示：

#生成试样
def GenSample(r,pos):
   for p in pos:
      body = GJKPolyhedron([],[1.e-2,1.e-2,1.e-2],0.05*1e-2,p_mat,False)
      body.state.pos=p
      O.bodies.append(body)
      O.bodies[-1].shape.color=(rand.random(), rand.random(), rand.random())

# 生成栅格
pos = GridInitial(R,num_x=num_x,num_y=num_y,num_z=num_z) # get positions of all nodes
# 为每个节点生成颗粒
GenSample(R,pos)

图 8.11 没有设置随机初始方向的多面体颗粒的最终堆积状态 (未显示立方体盒子)

图 8.12 没有设置随机初始方向的圆锥体的最终堆积状态 (未显示立方体盒子)

对于圆锥体来说函数 GenSample 需做如下改动:

#生成试样
def GenSample(r,pos):
   for p in pos:
      body = GJKCone(0.005,0.01,0.05*1e-2,p_mat, False)
      body.state.pos=p
      O.bodies.append(body)
      O.bodies[-1].shape.color=(rand.random(), rand.random(), rand.random())

图 8.12 给出了固定初始朝向的圆锥体试样 (通过将变量 rotate 设置为 False 来取消随机的初始朝向生成)。 可以看到所有颗粒的初始状态都是完全相同的。

用户可以尝试将变量 rotate 设置为 True 来测试下随机初始朝向的颗粒试样:

body = GJKCone(0.005,0.01,0.05*1e-2,p_mat,True)

图 8.13 设置了随机初始方向的圆锥体的最终堆积状态 (未显示立方体盒子)

如 图 8.13 所示，自由落体过程破坏了颗粒非常一致的初始朝向分布：

图 8.14 没有设置随机初始方向的圆柱体的最终堆积状态 (未显示立方体盒子)

圆柱的生成函数 GenSample 改动如下:

def GenSample(r,pos):
   for p in pos:
      body = GJKCylinder(0.005,0.01,0.05*1e-2,p_mat,False)
      body.state.pos=p
      O.bodies.append(body)
      O.bodies[-1].shape.color=(rand.random(), rand.random(), rand.random())

图 8.15 设置了随机初始方向的圆柱体的最终堆积状态 (未显示立方体盒子)

不出意料，对于初始朝向一致的圆柱颗粒，即使是在自由落体后，它们的朝向仍然不变 (图 8.14)。但是当初始朝向是随机生成时，自由落体结束后的圆柱体朝向变得不规律起来 (图 8.15)。

body = GJKCylinder(0.005,0.01,0.05*1e-2,p_mat,True)

图 8.16 没有设置随机初始方向的立方体的最终堆积状态 (未显示立方体盒子)

图 8.17 设置了随机初始方向的立方体的最终堆积状态 (未显示立方体盒子)

立方体的函数如下:

#生成试样
def GenSample(r,pos):
   for p in pos:
      body = GJKCuboid([0.005,0.005,0.005],0.05*1e-2,p_mat,False)
      body.state.pos=p
      O.bodies.append(body)
      O.bodies[-1].shape.color=(rand.random(), rand.random(), rand.random())

立方体的初始一致或随机生成朝向的自由落体过程与圆柱体非常相似。在初始朝向一致时，如 图 8.16 自由落体结束后朝向并未发生显著变化。 但是随机初始朝向的试样，在自由落体过程前后的朝向分布有巨大的差异。 True as follows:

body = GJKCuboid([0.005,0.005,0.005],0.05*1e-2,p_mat,True)

另外，考虑到立方体颗粒间的间隙较小，所以它们的随机朝向在自由落体前后的分布差异较圆柱体不明显。