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
|
# 球体在重力作用下自由落体到地面
from sudodem import utils # module utils has some auxiliary functions
#1. 定义球体的材料
mat = RolFrictMat(label="mat1",Kn=1e8,Ks=7e7,frictionAngle=math.atan(0.0),density=2650) # material 1 for particles
wallmat1 = RolFrictMat(label="wallmat1",Kn=1e8,Ks=0,frictionAngle=0.) # material 2 for walls
# 将新定义的材料添加到模拟的材料数值中
# O 可被认为是当前模拟的核心类
# 可以看到材料是 O 的一个属性, 所以我们可以用点操作来获取到对应的材料对象。
O.materials.append(mat) # materials is a list, i.e., a container for all materials involving in the simulations
O.materials.append(wallmat1)
# 在将材料添加到材料数组后我们可以通过材料的标签来获取其信息。参见下面墙体的材料定义。
# 创建球体
sp = sphere((0,0,10.0),1.0, material = mat)
O.bodies.append(sp)
# 用无限大墙体来模拟地面
ground = utils.wall(0,axis=2,sense=1, material = 'wallmat1') # utils.wall defines a wall with normal parallel to one axis of the global Cartesian coordinate system. The first argument (here = 0) specifies the location of the wall, and with the keyword axis (here = 2, axis has a value of 0, 1, 2 for x, y, z axies respectively) we know that the normal of the wall is along z axis with centering at z = 0 (0, specified by the first argument); the second keyword sense (here = 1, sense can be -1, 0, 1 for negative, both, and positive sides respectively) specifies that the positive side (i.e., the normal points to the positive direction of the axis) of the wall is activated, at which we will compute the interaction of a particle and this wall; the last keyword materials specifies a material to this wall, and we assign a string (i.e., the label of wallmat1) to the keyword, and directly assigning the object of the material (i.e., wallmat1 returned by RolFrictMat) is also acceptable.
O.bodies.append(ground) # add the ground to the body container
# 使用PyRunner来周期性地执行函数。
def record_data():
#open a file
fout = open('data.dat','a') #create/open a file named 'data.dat' in append mode
#we want to record the position of the sphere
p = O.bodies[0] # the sphere has been appended into the list O.bodies, and we can use the corresponding index to access it. We have added two bodies in total to O.bodies, and the sphere is the first one, i.e., the index is 0.
# p is the object of the sphere
# as a body object, p has several properties, e.g., state, material, which can be listed by p.dict() in the commandline.
pos = p.state.pos # get the position of the sphere
print>>fout, O.iter, pos[0], pos[1], pos[2] # we print the iteration number and the position (x,y,z) into the file fout.
# then, close the file and release resource
fout.close()
# 创建牛顿引擎来使得球体的运动遵循牛顿定律
newton=NewtonIntegrator(damping = 0.1,gravity=(0.,0.0,-9.8),label="newton")
# 设置局部阻尼为 0.1, 同时设置重力加速度为 9.8 m/s2，方向为z轴朝下。
# 引擎数组是整个模拟的最核心的部分
# 在每个时步内，引擎依次运行来完成 DEM 一个完整的计算循环
O.engines=[
ForceResetter(), # first, we need to reset the force container to store upcoming contact forces
# next, we conduct the broad phase of contact detection by comparing axis-aligned bounding boxs (AABBs) to rule out those pairs that are definitely not contacting.
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Wall_Aabb()],verletDist=0.2*0.01),
# then, we execute the narrow phase of contact detection
InteractionLoop( # we will loop all potentially contacting pairs of particles
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Wall_Sphere_ScGeom()], # for different particle shapes, we will call the corresponding functions to compute the contact geometry, for example, Ig2_Sphere_Sphere_ScGeom() is used to compute the contact geometry of a sphere-sphere pair.
[Ip2_RolFrictMat_RolFrictMat_RolFrictPhys()], # we will compute the physical properties of contact, e.g., contact stiffness and coefficient of friction, according to the physical properties of contacting particles.
[RollingResistanceLaw(use_rolling_resistance=False)] # Next, we compute contact forces in terms of the contact law with the info (contact geometric and physical properties) computed at last two steps.
),
newton, # finally, we compute acceleration and velocity of each particle and update their positions.
#at each time step, we have a PyRunner function help to hack into a DEM cycle and do whatever you want, e.g., changing particles' states, and saving some data, or just stopping the program after a certain running time.
PyRunner(command='record_data()',iterPeriod=10000,label='record',dead = False)
]
# 设置时步
O.dt = 1e-5
# 清空 data.dat 并写入文件头部信息
fout = open('data.dat','w') # we open the file in a write mode so that the file will be clean
print>>fout, 'iter, x, y, z' # write a file head
fout.close() # close the file
# 用户可通过如下两种方式来运行该模拟
# 1. 执行命令
# O.run()
# 或者设置时步来运行
O.run(1600000)
# 2. 在 GUI 界面点击运行按钮来运行
|
