tshortening.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
---
tshortening.py (5354B)
---
1 #!/usr/bin/env python
2 import sphere
3 import numpy
4
5 cube = sphere.sim('cube-init')
6 cube.readlast()
7 cube.adjustUpperWall(z_adjust=1.0)
8
9 # Fill out grid with cubic packages
10 grid = numpy.array((
11 [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
12 [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
13 [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
14 [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
15 [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
16 [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
17 [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
18 [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2],
19 [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1],
20 [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2],
21 [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
22 [ 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2],
23 [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
24 [ 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2],
25 [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
26 [ 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2],
27 [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
28 [ 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2],
29 [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
30 ))
31
32 # World dimensions and cube grid
33 nx = 1 # horizontal (thickness) cubes
34 ny = grid.shape[1] # horizontal cubes
35 nz = grid.shape[0] # vertical cubes
36 dx = cube.L[0]
37 dy = cube.L[1]
38 dz = cube.L[2]
39 Lx = dx*nx
40 Ly = dy*ny
41 Lz = dz*nz
42
43 sim = sphere.sim('shortening-relaxation', nw=0)
44
45 # insert particles into each cube in 90 degree CCW rotated coordinate system
46 # around y
47 for z in range(nz):
48 for y in range(ny):
49 for x in range(nx):
50
51 if (grid[z,y] == 0):
52 continue # skip to next iteration
53
54 for i in range(cube.np):
55 # x=x, y=Ly-z, z=y
56 pos = [ cube.x[i,0] + x*dx,
57 Ly - ((dz - cube.x[i,2]) + z*dz),
58 cube.x[i,1] + y*dy ]
59 sim.addParticle(pos, radius=cube.radius[i], color=grid[z,y])
60
61 # move to x=0
62 min_x = numpy.min(sim.x[:,0] - sim.radius[:])
63 sim.x[:,0] = sim.x[:,0] - min_x
64
65 # move to y=0
66 min_y = numpy.min(sim.x[:,1] - sim.radius[:])
67 sim.x[:,1] = sim.x[:,1] - min_y
68
69 # move to z=0
70 min_z = numpy.min(sim.x[:,2] - sim.radius[:])
71 sim.x[:,2] = sim.x[:,2] - min_z
72
73 #sim.defineWorldBoundaries(L=[Lx, Lz*3, Ly])
74 sim.defineWorldBoundaries(L=[numpy.max(sim.x[:,0] + sim.radius[:]), Lz*3, Ly])
75 sim.k_t[0] = 2.0/3.0*sim.k_n[0]
76
77 #sim.cleanup()
78 sim.writeVTK()
79 print(sim.np)
80
81
82 ## Relaxation
83
84 # Add gravitational acceleration
85 # Flip geometry so the upper wall pushes downwards
86 sim.g[0] = 0
87 sim.g[1] = -9.81
88 sim.g[2] = 0
89
90 sim.setDampingNormal(5.0e1)
91 sim.setDampingTangential(1.0e1)
92
93 #sim.periodicBoundariesX()
94 sim.normalBoundariesXY()
95 sim.uniaxialStrainRate(wvel = 0.0)
96
97 # Set duration of simulation, automatically determine timestep, etc.
98 sim.initTemporal(total=3.0, file_dt = 0.01)
99 sim.zeroKinematics()
100
101 sim.run(dry=True)
102 sim.run()
103 sim.writeVTKall()
104
105
106 ## Shortening
107 sim = sphere.sim('shortening-relaxation', nw=1)
108 sim.readlast()
109 sim.sid = 'shortening'
110 sim.cleanup()
111 sim.initTemporal(current=0.0, total=20.0, file_dt = 0.01, epsilon=0.07)
112
113 # set colors again
114 y_min = numpy.min(sim.x[:,1])
115 y_max = numpy.max(sim.x[:,1])
116 z_min = numpy.min(sim.x[:,2])
117 z_max = numpy.max(sim.x[:,2])
118 color_ny = 6
119 color_nz = int((z_max - z_min)/(y_max - y_min)*color_ny)
120 color_dy = y_max/color_ny
121 color_dz = z_max/color_nz
122 color_y = numpy.arange(0.0, y_max, ny)
123 color_z = numpy.arange(0.0, z_max, nz)
124
125 # 1 or 2 in horizontal layers
126 #for i in range(ny-1):
127 #I = numpy.nonzero((sim.x[:,1] >= color_y[i]) & (sim.x[:,1] <= color_y[i+1]))
128 #sim.color[I] = i%2 + 1
129
130 # 1 or 3 in checkerboard
131 for i in range(sim.np):
132 iy = numpy.floor((sim.x[i,1] - y_min)/(y_max/color_ny))
133 iz = numpy.floor((sim.x[i,2] - z_min)/(z_max/color_nz))
134 sim.color[i] = (-1)**iy + (-1)**iz + 1
135
136 # fix lowest plane of particles
137 I = numpy.nonzero(sim.x[:,1] < 1.5*numpy.mean(sim.radius))
138 sim.fixvel[I] = -1
139 sim.color[I] = 0
140 sim.x[I,1] = 0.0 # move into center into lower wall to avoid stuck particles
141
142 # fix left-most plane of particles
143 I = numpy.nonzero(sim.x[:,2] < 1.5*numpy.mean(sim.radius))
144 sim.fixvel[I] = -1
145 sim.color[I] = 0
146
147 # fix right-most plane of particles
148 I = numpy.nonzero((sim.x[:,2] > z_max - 1.5*numpy.mean(sim.radius)) &
149 (sim.x[:,1] > 10.0*numpy.mean(sim.radius)))
150 sim.fixvel[I] = -1
151 sim.color[I] = 0
152
153 #sim.normalBoundariesXY()
154 #sim.periodicBoundariesX()
155 sim.zeroKinematics()
156
157 # Wall parameters
158 sim.mu_ws[0] = 0.5
159 sim.mu_wd[0] = 0.5
160 #sim.gamma_wn[0] = 1.0e2
161 #sim.gamma_wt[0] = 1.0e2
162 sim.gamma_wn[0] = 0.0
163 sim.gamma_wt[0] = 0.0
164
165 # Particle parameters
166 sim.mu_s[0] = 0.5
167 sim.mu_d[0] = 0.5
168 sim.gamma_n[0] = 1000.0
169 sim.gamma_t[0] = 1000.0
170
171 # push down upper wall
172 compressional_strain = 0.5
173 wall_velocity = -compressional_strain*(z_max - z_min)/sim.time_total[0]
174 sim.uniaxialStrainRate(wvel = wall_velocity)
175 sim.vel[I,2] = wall_velocity
176
177 sim.run(dry=True)
178 sim.run()
179 sim.writeVTKall()