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
92
93
94
|
##################################################
from sudodem import utils
from sudodem._superellipse_utils import *
import numpy as np
import math
import random
random.seed(11245642)
####################################################
#########设置
####################################################
isSphere = False
r = 0.5 #颗粒大小
num_s=100
num_t=8000
trials=0
boxsize=[50.0,10.0]
#############################
mat =SuperellipseMat(label="mat1",Kn=1e8,Ks=7e7,frictionAngle=math.atan(0.225),density=2650)
O.materials.append(mat)
O.bodies.append(utils.wall(-0.5*boxsize[1],axis=1,sense=1, material = 'mat1'))#ground
O.bodies.append(utils.wall(-0.5*boxsize[0],axis=0,sense=1, material = 'mat1'))#left
O.bodies.append(utils.wall(0.5*boxsize[0],axis=0,sense=-1, material = 'mat1'))#right
###
def gettop():
ymax=0
for b in O.bodies:
if(isinstance(b.shape,Superellipse)):
if(b.state.pos[1]>=ymax):
ymax=b.state.pos[1]
return ymax
def Genparticles(r,mat,boxsize,num_s,num_t,bottom):
global trials
gap=r
#print(gap)
num=0
coor=[]
width=(boxsize[0]-2.*gap)/2.
height=(boxsize[1]-2.*gap)
#iteration=0
while num<num_s and trials<num_t:
isOK=True
pos=[0]*2
pos[0]=random.uniform(-width,width)
pos[1]=random.uniform(bottom+gap,bottom+gap+height)
for i in range(0,len(coor)):
distance=sum([((coor[i][j]-pos[j])**2.) for j in range(0,2)])
if(distance<(4.*gap*gap)):
isOK=False
break
if(isOK==True):
coor.append(pos)
num+=1
trials+=1
#print(num)
return coor
def Addlayer(r,mat,boxsize,num_s,num_t):
bottom=gettop()
coor=Genparticles(r,mat,boxsize,num_s,num_t,bottom)
for b in coor:
rr = r*random.uniform(0.5,1.0)
bb = NewSuperellipse(1.5*r,r,1.0,mat,True,isSphere)
bb.state.pos=b
bb.state.vel=[0,-1]
O.bodies.append(bb)
if len(O.bodies) - 3 > 2000:
O.engines[-1].dead=True
O.dt = 1e-3
newton=NewtonIntegrator(damping = 0.1,gravity=(0.,-10.0),label="newton",isSuperellipse=True)
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Superellipse_Aabb(),Bo1_Wall_Aabb()],verletDist=0.1),
InteractionLoop(
[Ig2_Wall_Superellipse_SuperellipseGeom(), Ig2_Superellipse_Superellipse_SuperellipseGeom()],
[Ip2_SuperellipseMat_SuperellipseMat_SuperellipsePhys()], # 接触检测引擎
[SuperellipseLaw()] # 接触力本构模型
),
newton,
PyRunner(command='Addlayer(r,mat,boxsize,num_s,num_t)',virtPeriod=0.1,label='check',dead = False)
]
|
{#figpackingellipse width=“10cm”}