6.2. 颗粒柱塌落

本算例展示如何使用CuDEM模拟颗粒柱在重力作用下的塌落。如 图 6.2 所示，参照上一个例子 颗粒堆积 ，首先生成颗粒柱，把容器扩大，随后颗粒柱在重力作用下塌落。需要指出的是，本算例将生成的颗粒柱并未在重力作用下先平衡（即上一例子的最终堆积状态）。如你感兴趣，可以采用堆积后的颗粒柱进行塌落模拟。

../_images/columncollapsesetup.png

图 6.2 160万颗粒的塌落柱模拟

6.2.1. 导入基础包

导入pycudem和pySudoMath:

1import sys,os,math,random
2import numpy as np
3from dem.core import *
4from dem.cudem import *
5from pySudoMath import *

6.2.2. 设置GUI或命令行运行

为了让脚本既能在GUI中运行，又能直接通过命令行运行，我们做如下设置：

1gui = True
2sim = None
3if gui:
4    app.core.setAnalysisType(AnalysisType="CUDEM")
5    sim=app.core.sim
6else:
7    sim = DEMSim()

6.2.3. 定义一些辅助函数

  • 定义生成栅格的函数，我们将在栅格点的位置生成颗粒：

 1#define a lattice grid
 2#R: distance between two neighboring nodes
 3#num_x: number of nodes along x axis
 4#num_y: number of nodes along y axis
 5#num_z: number of nodes along z axis
 6#return a list of the postions of all nodes
 7def GridInitial(R,num_x=10,num_y=10,num_z=20, shift=[0,0,0]):#this part is slow, but just for test
 8    pos = list()
 9    for k in range(num_z):
10        for j in range(num_y):
11            for i in range(num_x):
12                x = i*R*2.0*1.2+R+shift[0]+random.uniform(0, 0.1*R)
13                y = j*R*2.0*1.2+R+shift[1]+random.uniform(0, 0.1*R)
14                z = k*R*2.0*1.2+R+shift[2]
15                pos.append(Vector3f(x,y,z))
16    return pos

  • 定义生成墙的函数

 1def make_wall(R,num_x=1,num_y=1,num_z=1, ll_corner=Vector3f(0,0,0)):#this part is slow, but just for test
 2    pos = list()
 3    for k in range(num_z):
 4        for j in range(num_y):
 5            for i in range(num_x):
 6                x = i*R*2.0+R
 7                y = j*R*2.0+R
 8                z = k*R*2.0+R
 9                pos.append(Vector3f(x,y,z)+ll_corner)
10    return pos

6.2.4. 设置模型大小

设置颗粒AABB包围盒的放大系数 skin，设置是否启用光追核心(RT)加速RTflag，设置颗粒大小（单位m）及模拟盒子的大小（颗粒大小的倍数）：

1#skin size
2skin = 1.2
3RTflag = False#True
4#modelScale 1 : 1.6M particles
5modelScale = 16
6psize = 0.1
7boxsize_x = 71
8boxsize_y = 10*modelScale
9boxsize_z = 180

6.2.5. 设置墙

  • 设置前后左右及底部墙：

1#walls
2sim.walls.addWall(Vector3f(0,0,0), 1, 2) #bottom
3sim.walls.addWall(Vector3f(0,0,0), 1, 0) #left
4sim.walls.addWall(Vector3f(boxsize_x*psize*2*1.2+psize*2,0,0), -1, 0) #right
5sim.walls.addWall(Vector3f(0,0,0), 1, 1) #front
6sim.walls.addWall(Vector3f(0,boxsize_y*psize*2*1.2+psize*2,0), -1, 1) #back

  • 可以通过索引获得生成的墙，比如：

1wall_right = sim.walls[2]

  • 设置墙的材料，采用线弹性模型

1wmat = Material()
2wmat.kn = 1e7 #法向刚度
3wmat.ks = 1e7 #切向刚度
4wmat.friction = 1.0 #摩擦系数
5wmat.vdamping = 0.5 #粘滞阻尼
6sim.walls.material = wmat #把材料赋值到所有墙

6.2.6. 设置模拟域

设置模拟域的AABB大小以及模拟域中分域大小（采用基于cell的算法需要此设置）：

1#domain size should be large enough
2sim.domain.setAABB(psize*boxsize_x*2*1.2*2+psize*2,psize*boxsize_y*2*1.2+psize*2,psize*boxsize_z*2)
3sim.domain.setCellSize(psize*2.0*skin)

6.2.7. 设置默认材料

颗粒采用默认材料，参考墙的材料：

1sim.material.kn = 1e5
2sim.material.ks = 1e5
3sim.material.friction = 0.5
4sim.material.density = 2000.0
5sim.material.vdamping = 0.5

6.2.8. 设置运动求解器

设置阻尼，重力加速度及 alpha 。其中， alpha 将会影响邻居列表的更新频率。

1sim.integrator.damping = 0.1
2sim.integrator.gravity = Vector3f(0.0, 0.0, -10.0)
3sim.integrator.alpha = skin

6.2.9. 设置邻居列表属性

颗粒运动超过 delta 倍颗粒粒径后更新邻居列表； 光追接触检测时放大的AABB大小为 ORTextend 。

1sim.nlist.delta = skin
2sim.nlist.ORTextend = 2.0*psize*skin

6.2.10. 设置内存动态优化

设置颗粒重新内存从新排列频率 reorderInterval （以时步计）；设置是否开启光追接触计算 useORT ； 设置计算时间步长 dt 。

1sim.reorderInterval = 500 #2000
2sim.useORT = RTflag
3sim.dt = 1e-3

6.2.11. 生成颗粒

  • 生成颗粒位置和大小，本例子采用单一粒径颗粒：

 1Num_x = 71
 2Num_y = 10*modelScale
 3Num_z = 141#120
 4pos = GridInitial(psize,num_x=Num_x,num_y=Num_y,num_z=Num_z,shift = [1.0*psize,1.0*psize,1.0*psize])
 5
 6#rs = np.random.rand(Num*Num*Num)*0.5+0.5
 7#rs = np.ones(Num*Num*Num_z+len(wall_pos))*psize
 8rs_p = np.ones(Num_x*Num_y*Num_z)*psize#*random.uniform(0.5,1.0)
 9# rs_w = np.ones(len(wall_pos))*psize
10rs = rs_p #np.hstack((rs_w,rs_p))

  • 根据颗粒位置和半径生成颗粒：

1sim.createParticles(pos,rs)

6.2.12. 设置钩子程序

  • 设置辅助函数 HookFunc ，该辅助函数用于记录运行时间，可以在模拟计算中自动调用。

1import time
2elapsedTime = []
3def HookFunc(elapsedTime):
4    t = time.time()
5    if len(elapsedTime) > 0:
6        elapsedTime.append(t-elapsedTime[0])
7    else:
8        elapsedTime.append(t)

  • 设置钩子程序，并加到模拟中。

1pyhook1 = PyHook()
2# 指定该钩子程序需要执行的Python函数
3pyhook1.command = 'HookFunc(elapsedTime)'
4# 设置每隔800步执行一次，总共执行100次
5pyhook1.reset(iterPeriod = 800, totalRuns = 100)
6# 让该钩子程序处于执行状态
7pyhook1.dead = False
8# 添加该钩子到模拟中
9sim.hooks = [pyhook1]

6.2.13. 运行程序

  • GUI中命令行

1app.run(10000)

  • 终端命令行

1sudosim -i -f packing.py

然后在Python解释器中运行：

1sim.RunPthread(10000)