8.7. 超椭圆堆积模拟

该示例给出了超椭圆的堆积模拟 (超椭球的二维形状)

8.7.1. 超椭圆

超椭圆局部坐标系下的形状函数如下式：

\[\big|\frac{x}{r_x}\big|^{\frac{2}{\epsilon}} + \big|\frac{y}{r_y}\big|^{\frac{2}{\epsilon}} = 1\]

其中 \(r_x\) 和 \(r_y\) 分别表示x和y轴的半长轴长度; \(\epsilon \in (0, 2)\) 控制了颗粒表面的棱角度。

8.7.2. 建模脚本

我们这里给出了超椭圆在重力作用下自由落体堆积的模拟脚本:

  1. 导入基础包
    1from sudodem import utils
2from sudodem._superellipse_utils import *
3import numpy as np
4import math
5import random
6random.seed(11245642)
  2. 模拟参数设置
    1isSphere = False
2r = 0.5 #颗粒大小
3num_s=100
4num_t=8000
5trials=0
6boxsize=[50.0,10.0]
  3. 生成墙/容器
    1mat =SuperellipseMat(label="mat1",Kn=1e8,Ks=7e7,frictionAngle=math.atan(0.225),density=2650)
2O.materials.append(mat)
3O.bodies.append(utils.wall(-0.5*boxsize[1],axis=1,sense=1, material = 'mat1'))#ground
4O.bodies.append(utils.wall(-0.5*boxsize[0],axis=0,sense=1, material = 'mat1'))#left
5O.bodies.append(utils.wall(0.5*boxsize[0],axis=0,sense=-1, material = 'mat1'))#right
  4. 定义辅助函数
     1def gettop():
 2    ymax=0
 3    for b in O.bodies:
 4        if(isinstance(b.shape,Superellipse)):
 5            if(b.state.pos[1]>=ymax):
 6                ymax=b.state.pos[1]
 7    return ymax
 8
 9
10def Genparticles(r,mat,boxsize,num_s,num_t,bottom):
11    global trials
12    gap=r
13    #print(gap)
14    num=0
15    coor=[]
16    width=(boxsize[0]-2.*gap)/2.
17    height=(boxsize[1]-2.*gap)
18    #iteration=0
19
20    while num<num_s and trials<num_t:
21        isOK=True
22        pos=[0]*2
23        pos[0]=random.uniform(-width,width)
24        pos[1]=random.uniform(bottom+gap,bottom+gap+height)
25
26        for i in range(0,len(coor)):
27            distance=sum([((coor[i][j]-pos[j])**2.) for j in range(0,2)])
28            if(distance<(4.*gap*gap)):
29
30                isOK=False
31                break
32
33        if(isOK==True):
34            coor.append(pos)
35            num+=1
36            trials+=1
37            #print(num)
38    return coor
39
40
41def Addlayer(r,mat,boxsize,num_s,num_t):
42    bottom=gettop()
43    coor=Genparticles(r,mat,boxsize,num_s,num_t,bottom)
44    for b in coor:
45        rr = r*random.uniform(0.5,1.0)
46        bb = NewSuperellipse(1.5*r,r,1.0,mat,True,isSphere)
47        bb.state.pos=b
48        bb.state.vel=[0,-1]
49        O.bodies.append(bb)
50    if len(O.bodies) - 3 > 2000:
51        O.engines[-1].dead=True
  5. 模拟引擎设置
     1newton=NewtonIntegrator(damping = 0.1,gravity=(0.,-10.0),label="newton",isSuperellipse=True)
 2
 3O.engines=[
 4ForceResetter(),
 5InsertionSortCollider([Bo1_Superellipse_Aabb(),Bo1_Wall_Aabb()],verletDist=0.1),
 6InteractionLoop(
 7    [Ig2_Wall_Superellipse_SuperellipseGeom(), Ig2_Superellipse_Superellipse_SuperellipseGeom()],
 8    [Ip2_SuperellipseMat_SuperellipseMat_SuperellipsePhys()], # 接触检测引擎
 9    [SuperellipseLaw()]   # 接触力本构模型
10),
11newton,
12PyRunner(command='Addlayer(r,mat,boxsize,num_s,num_t)',virtPeriod=0.1,label='check',dead = False)
13]
14
15O.dt = 1e-3

如上脚本需使用 sudodem2d 运行，模拟结果如 图 8.20。

../_images/packingellipse.png

图 8.20 超椭圆在重力作用下的堆积状态