tchannel.py - sphere - GPU-based 3D discrete element method algorithm with optional fluid coupling
(HTM) git clone git://src.adamsgaard.dk/sphere
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tchannel.py (5432B)
---
1 #!/usr/bin/env python
2 import sphere
3 import numpy
4
5 relaxation = True
6 consolidation = False
7 water = False
8
9 id_prefix = 'channel3'
10 N = 10e3
11
12 cube = sphere.sim('cube-init')
13 cube.readlast()
14 cube.adjustUpperWall(z_adjust=1.0)
15
16 # Fill out grid with cubic packages
17 grid = numpy.array((
18 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
19 [1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
20 [1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
21 [1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
22 [1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
23 [1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
24 [1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
25 [1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
26 [1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
27 [1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
28 [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
29 [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
30 [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
31 [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
32 [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
33 [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
34 [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
35 [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
36 [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
37 [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
38 [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
39 [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
40 [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
41 [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
42 [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
43 [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
44 [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
45 [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
46 [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
47 [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
48 ))
49
50 # World dimensions and cube grid
51 nx = grid.shape[1] # horizontal cubes
52 ny = 1 # horizontal (thickness) cubes
53 nz = grid.shape[0] # vertical cubes
54 dx = cube.L[0]
55 dy = cube.L[1]
56 dz = cube.L[2]
57 Lx = dx*nx
58 Ly = dy*ny
59 Lz = dz*nz
60
61 sim = sphere.sim(id_prefix + '-relaxation', nw=0)
62
63 # insert particles into each cube
64 for z in range(nz):
65 for y in range(ny):
66 for x in range(nx):
67
68 if (grid[z, x] == 0):
69 continue # skip to next iteration
70
71 for i in range(cube.np):
72 pos = [cube.x[i, 0] + x*dx,
73 cube.x[i, 1] + y*dy,
74 cube.x[i, 2] + (nz-z)*dz]
75
76 sim.addParticle(pos, radius=cube.radius[i], color=grid[z, x])
77
78 # move to x=0
79 min_x = numpy.min(sim.x[:, 0] - sim.radius[:])
80 sim.x[:, 0] = sim.x[:, 0] - min_x
81
82 # move to y=0
83 min_y = numpy.min(sim.x[:, 1] - sim.radius[:])
84 sim.x[:, 1] = sim.x[:, 1] - min_y
85
86 # move to z=0
87 min_z = numpy.min(sim.x[:, 2] - sim.radius[:])
88 sim.x[:, 2] = sim.x[:, 2] - min_z
89
90 sim.defineWorldBoundaries(L=[numpy.max(sim.x[:, 0] + sim.radius[:]),
91 numpy.max(sim.x[:, 1] + sim.radius[:]),
92 numpy.max(sim.x[:, 2] + sim.radius[:])*10.2])
93 #numpy.max(sim.x[:, 2] + sim.radius[:])*1.2])
94 sim.k_t[0] = 2.0/3.0*sim.k_n[0]
95
96 # sim.cleanup()
97 sim.writeVTK()
98 print("Number of particles: " + str(sim.np))
99
100
101 # Relaxation
102
103 # Add gravitational acceleration
104 sim.g[0] = 0.0
105 sim.g[1] = 0.0
106 sim.g[2] = -9.81
107
108 sim.normalBoundariesXY()
109 # sim.consolidate(normal_stress=0.0)
110
111 # assign automatic colors, overwriting values from grid array
112 sim.checkerboardColors(nx=grid.shape[1], ny=2, nz=grid.shape[0]/4)
113
114 sim.contactmodel[0] = 2
115 sim.mu_s[0] = 0.5
116 sim.mu_d[0] = 0.5
117
118 # Set duration of simulation, automatically determine timestep, etc.
119 sim.initTemporal(total=3.0, file_dt=0.01, epsilon=0.07)
120 sim.time_dt[0] = 1.0e-20
121 sim.time_file_dt = sim.time_dt
122 sim.time_total = sim.time_file_dt*5.
123 sim.zeroKinematics()
124
125 if relaxation:
126 sim.run(dry=True)
127 sim.run()
128 sim.writeVTKall()
129
130 exit()
131
132 # Consolidation under constant normal stress
133 if consolidation:
134 sim.readlast()
135 sim.id(id_prefix + '-' + str(int(N/1000.)) + 'kPa')
136 sim.cleanup()
137 sim.initTemporal(current=0.0, total=10.0, file_dt=0.01, epsilon=0.07)
138
139 # fix lowest plane of particles
140 I = numpy.nonzero(sim.x[:, 2] < 1.5*numpy.mean(sim.radius))
141 sim.fixvel[I] = -1
142 sim.color[I] = 0
143
144 sim.zeroKinematics()
145
146 # Wall parameters
147 sim.mu_ws[0] = 0.5
148 sim.mu_wd[0] = 0.5
149 sim.gamma_wn[0] = 1.0e2
150 sim.gamma_wt[0] = 1.0e2
151 # sim.gamma_wn[0] = 0.0
152 # sim.gamma_wt[0] = 0.0
153
154 # Particle parameters
155 sim.mu_s[0] = 0.5
156 sim.mu_d[0] = 0.5
157 sim.gamma_n[0] = 0.0
158 sim.gamma_t[0] = 0.0
159
160 # apply effective normal stress from upper wall
161 sim.consolidate(normal_stress=N)
162
163 sim.run(dry=True)
164 #sim.run()
165 #sim.writeVTKall()
166
167 ## Add water
168 #if water:
169 #sim.readlast()
170 #sim.id(id_prefix + '-wet')
171 #sim.wet()
172 #sim.initTemporal(total=3.0, file_dt=0.01, epsilon=0.07)
173 #
174 #sim.run(dry=True)
175 #sim.run()
176 #sim.writeVTKall()