Introduction
Introduction Statistics Contact Development Disclaimer Help
tshortening.py - sphere - GPU-based 3D discrete element method algorithm with o…
git clone git://src.adamsgaard.dk/sphere
Log
Files
Refs
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 sy…
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…
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…
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…
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 parti…
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()
You are viewing proxied material from mx1.adamsgaard.dk. The copyright of proxied material belongs to its original authors. Any comments or complaints in relation to proxied material should be directed to the original authors of the content concerned. Please see the disclaimer for more details.