示例 1: 超级椭球堆积模拟

超级椭球局部坐标下的形状函数可由下式给出 $$\label{eqsurface} \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$$ 其中 $r_x, r_y$ 和 $r_z$ 是超级椭球在x, y, z轴的半长主轴; and $\epsilon_i (i=1,2)$ 是刻画颗粒棱角度的形状参数. $\epsilon_i$ 在 0 到 2 的变化可勾勒出非常丰富的凸超级椭球形状。SudoDEM 生成超级椭球的两个函数可参见如下高亮描述 (see Sec. 5.1.2{reference-type=“ref” reference=“modsuperquadrics”}):

• NewSuperquadrics2($r_x$, $r_y$, $r_z$, $\epsilon_1$, $\epsilon_2$, mat, rotate, isSphere)

• NewSuperquadrics_rot2($r_x$, $r_y$, $r_z$, $\epsilon_1$, $\epsilon_2$, mat, $q_w$, $q_x$, $q_y$, $q_z$, isSphere)

这里我们给出了超级椭球在立方体盒子中重力堆积的模拟示例。

 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

###############################################
# 在栅格中生成超级椭球，随后在重力作用下堆积在立方体中
###############################################
# 导入一些必要的库
from sudodem import _superquadrics_utils

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

######## 定义形状参数 #################
isSphere=False

num_x = 5
num_y = 5
num_z = 20
R = 0.1

####### 定义一些辅助函数 #########
               
# 定义栅格
# 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

# 生成试样
def GenSample(r,pos):
   for p in pos:
		epsilon1 = rand.uniform(0.7,1.6)
		epsilon2 = rand.uniform(0.7,1.6)
		a = r
		b = r*rand.uniform(0.4,0.9)#aspect ratio
		c = r*rand.uniform(0.4,0.9)
		
		body = NewSuperquadrics2(a,b,c,epsilon1,epsilon2,p_mat,True,isSphere)
		body.state.pos=p
		O.bodies.append(body)

######### 模拟设置 ####################
# 颗粒材料设置
p_mat = SuperquadricsMat(label="mat1",Kn=1e5,Ks=7e4,frictionAngle=math.atan(0.3),density=2650,betan=0,betas=0)
# betan 和 betas 对应于接触的粘性阻尼, 默认值为0，表示无粘性阻尼.
# 侧壁墙的材料设置
wall_mat = SuperquadricsMat(label="wallmat",Kn=1e6,Ks=7e5,frictionAngle=0.0,betan=0,betas=0)
# 底部墙体的材料设置
wallmat_b = SuperquadricsMat(label="wallmat",Kn=1e6,Ks=7e5,frictionAngle=math.atan(1),betan=0,betas=0)
# 将定义的材料添加到模拟材料数组中
O.materials.append(p_mat)
O.materials.append(wall_mat)
O.materials.append(wallmat_b)

# 创建墙体对象并将其添加到模拟的形状列表中
O.bodies.append(utils.wall(-R,axis=0,sense=1, material = wall_mat))#left wall along x axis
O.bodies.append(utils.wall(2.0*R*num_x-R,axis=0,sense=-1, material = wall_mat))#right wall along x axis
O.bodies.append(utils.wall(-R,axis=1,sense=1, material = wall_mat))#front wall along y axis
O.bodies.append(utils.wall(2.0*R*num_y-R,axis=1,sense=-1, material = wall_mat))#back wall along y axis
O.bodies.append(utils.wall(-R,axis=2,sense=1, material =wallmat_b))#bottom wall along z axis

# 创建栅格
pos = GridInitial(R,num_x=num_x,num_y=num_y,num_z=num_z) # get positions of all nodes
# 在每个节点创建超级椭球颗粒
GenSample(R,pos)

# 创建引擎
# 定义牛顿引擎
newton=NewtonIntegrator(damping = 0.1,gravity=(0.,0.,-9.8),label="newton",isSuperquadrics=1) # isSuperquadrics: 1 for superquadrics

O.engines=[
   ForceResetter(), InsertionSortCollider([Bo1_Superquadrics_Aabb(),Bo1_Wall_Aabb()],verletDist=0.2*0.1),
   InteractionLoop(
      [Ig2_Wall_Superquadrics_SuperquadricsGeom(),    Ig2_Superquadrics_Superquadrics_SuperquadricsGeom()], 
      [Ip2_SuperquadricsMat_SuperquadricsMat_SuperquadricsPhys()], # collision "physics"
      [SuperquadricsLaw()]   # contact law
   ),
   newton
   
]
O.dt=5e-5

在GUI控制面板中我们可以通过点击 ‘show 3D’ 按钮来显示模拟的可视化场景。在显示模式下选择 ‘Gl1_Superquadrics’ 渲染器来控制颗粒是由线条还是由实体面来渲染（如图）， 3.3{reference-type=“ref” reference=“figguidisplay”}.

Figs. 3.4{reference-type=“ref” reference=“figexp11”} - 3.6{reference-type=“ref” reference=“figexp13”} 显示了在不同渲染模式下的超级椭球颗粒。

超级椭球颗粒的属性和类方法展示如下：

• 属性:

  • $r_x$, $r_y$, $r_z$: 浮点型, x, y, 和 z 轴的半长轴。

  • $\epsilon_1$, $\epsilon_2$: 浮点型, 材料参数。

  • isSphere: 布尔值, 颗粒是否是球体.

• 类方法:

  • getVolume(): 浮点型, 返回颗粒体积。

  • getrxyz(): 三阶向量, 返回颗粒的三个半长轴。

  • geteps(): 二阶向量, 返回颗粒的形状参数eps1和eps2。

如下代码给出了获取索引为100的超级椭球的体积:

volume = O.bodies[100].shape.getVolume()

获取对象所有属性的方法如下:

O.bodies[100].dict()

输出如下:

{'bound': <Aabb instance at 0x56072ec83da0>,
 'chain': -1,
 'clumpId': -1,
 'flags': 3,
 'groupMask': 1,
 'id': 100,
 'iterBorn': 0,
 'material': <SuperquadricsMat instance at 0x56072e6ed5f0>,
 'shape': <Superquadrics instance at 0x56072ec83b60>,
 'state': <State instance at 0x56072ec839e0>,
 'timeBorn': 0.0}

我们可以看到对象的属性是以字典的方式被打印出来的。需要提及的是尖括号内的是另一个对象的指针，也就是说我们可以通过dict()方法来进一步获取其对应的属性。例如我们要获取颗粒state的属性，可以通过如下代码：

O.bodies[100].state.dict()

输出如下:

{'angMom': Vector3(0,0,0),
 'angVel': Vector3(0,0,0),
 'blockedDOFs': 0,
 'densityScaling': 1.0,
 'inertia': Vector3(0.0045389863901528615,0.007287422614611933,
                           0.0067053718092146865),
 'isDamped': True,
 'mass': 3.225715855242557,
 'refOri': Quaternion((1,0,0),0),
 'refPos': Vector3(0,0,0),
 'se3': (Vector3(0,0.8,3),
  Quaternion((0.3970010520699211,0.5339678220536823,
                     -0.7465042060609055),2.1754719537556495)),
 'vel': Vector3(0,0,0)}

对于 Shape对象, 它的属性输出可能随着子类的不同而不同。

O.bodies[100].shape.dict()

输出如下:

{'color': Vector3(0.6407663308273844,0.07426042997942327,
                            0.05872324344642611),
 'eps': Vector2(1.3975848801318045,1.5376309088540607),
 'eps1': 1.3975848801318045,
 'eps2': 1.5376309088540607,
 'highlight': False,
 'isSphere': False,
 'rx': 0.1,
 'rxyz': Vector3(0.1,0.06469580193313924,0.07572423080016706),
 'ry': 0.06469580193313924,
 'rz': 0.07572423080016706,
 'wire': False}

这里我们给出了一个例子用来获取所有超级椭球的序号及其位置信息:

for b in O.bodies: # loop all bodies in the simulation
    # we have to check if the present body is a superellipsoid or others, e.g., a wall
    if isinstance(b.shape, Superquadrics): # if it is a superellipsoid, then to do
        pid = b.id # get the id of particle b
        pos = b.state.pos # get the position of particle b
        print pid, pos