tsegregation.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
---
tsegregation.py (5552B)
---
1 #!/usr/bin/env python
2
3 # Import sphere functionality
4 import sphere
5
6 ### EXPERIMENT SETUP ###
7 initialization = True
8 consolidation = True
9 shearing = True
10 rendering = True
11 plots = True
12
13 # Number of particles
14 np = 2e4
15
16 # Common simulation id
17 sim_id = 'segregation'
18
19 # Deviatoric stress [Pa]
20 #devslist = [80e3, 10e3, 20e3, 40e3, 60e3, 120e3]
21 devslist = [120e3]
22 #devs = 0
23
24 ### INITIALIZATION ###
25
26 # New class
27 init = sphere.sim(np = np, nd = 3, nw = 0, sid = sim_id + '-init')
28
29 # Save radii
30 init.generateRadii(mean = 0.08)
31
32 # Use default params
33 init.defaultParams(gamma_n = 100.0, mu_s = 0.4, mu_d = 0.4)
34
35 init.periodicBoundariesXY()
36
37 # Initialize positions in random grid (also sets world size)
38 hcells = np**(1.0/3.0)*0.6
39 init.initRandomGridPos(gridnum = numpy.array([hcells, hcells, 1e9]))
40
41 # Choose the tangential contact model
42 # 1) Visco-frictional (somewhat incorrect, fast computations)
43 # 2) Elastic-viscous-frictional (more correct, slow computations in dense
44 # packings)
45 init.contactmodel[0] = 2
46
47 # Set horizontal boundaries as periodic
48 init.periodicBoundariesXY()
49
50 # Add gravity
51 init.g[2] = -9.81
52
53 # Decrease size of top particles
54 fraction = 0.80
55 z_min = numpy.min(init.x[:,2] - init.radius)
56 z_max = numpy.max(init.x[:,2] + init.radius)
57 I = numpy.nonzero(init.x[:,2] > (z_max - (z_max-z_min)*fraction))
58 init.radius[I] = init.radius[I] * 0.5
59
60 # Set duration of simulation
61 init.initTemporal(total = 10.0)
62
63 if (initialization == True):
64
65 # Run sphere
66 init.run(dry=True)
67 init.run()
68
69 if (plots == True):
70 # Make a graph of energies
71 init.visualize('energy', savefig=True, outformat='png')
72
73 if (rendering == True):
74 # Render images with raytracer
75 init.render(method = 'angvel', max_val = 0.3, verbose = False)
76
77
78 ### CONSOLIDATION ###
79
80 for devs in devslist:
81 # New class
82 cons = sphere.sim(np = init.np, nw = 1, sid = sim_id
83 + '-cons-devs{}'.format(devs))
84
85 # Read last output file of initialization step
86 lastf = sphere.status(sim_id + '-init')
87 cons.readbin('../output/' + sim_id
88 + '-init.output{:0=5}.bin'.format(lastf), verbose=False)
89
90 # Setup consolidation experiment
91 cons.consolidate(normal_stress = devs, periodic = init.periodic)
92
93 # Set duration of simulation
94 cons.initTemporal(total = 1.5)
95
96 if (consolidation == True):
97
98 # Run sphere
99 cons.run(dry=True) # show values, don't run
100 cons.run() # run
101
102 if (plots == True):
103 # Make a graph of energies
104 cons.visualize('energy', savefig=True, outformat='png')
105 cons.visualize('walls', savefig=True, outformat='png')
106
107 if (rendering == True):
108 # Render images with raytracer
109 cons.render(method = 'pres', max_val = 2.0*devs, verbose = False)
110
111
112 ### SHEARING ###
113
114 # New class
115 shear = sphere.sim(np = cons.np, nw = cons.nw,
116 sid = sim_id + '-shear-devs{}'.format(devs))
117
118 # Read last output file of initialization step
119 lastf = sphere.status(sim_id + '-cons-devs{}'.format(devs))
120 shear.readbin('../output/' + sim_id
121 + '-cons-devs{}.output{:0=5}.bin'.format(devs, lastf),
122 verbose = False)
123
124 # Setup shear experiment
125 #shear.shear(shear_strain_rate = 0.05, periodic = init.periodic)
126 shear_strain_rate = 0.05
127
128 ## Custom shear function
129 # Find lowest and heighest point
130 z_min = numpy.min(shear.x[:,2] - shear.radius)
131 z_max = numpy.max(shear.x[:,2] + shear.radius)
132
133 # the grid cell size is equal to the max. particle diameter
134 cellsize = shear.L[0] / shear.num[0]
135
136 # make grid one cell heigher to allow dilation
137 shear.num[2] += 1
138 shear.L[2] = shear.num[2] * cellsize
139
140 # zero kinematics
141 shear.zeroKinematics()
142
143 # set friction coefficients
144 shear.mu_s[0] = 0.5
145 shear.mu_d[0] = 0.5
146
147 # set the thickness of the horizons of fixed particles
148 #fixheight = 2*cellsize
149 #fixheight = cellsize
150
151 # Fix horizontal velocity to 0.0 of lowermost particles
152 d_max_below = numpy.max(shear.radius[numpy.nonzero(shear.x[:,2] < \
153 (z_max-z_min)*0.3)])*2.0
154 #I = numpy.nonzero(shear.x[:,2] < (z_min + fixheight))
155 I = numpy.nonzero(shear.x[:,2] < (z_min + d_max_below))
156 shear.fixvel[I] = 1
157 shear.angvel[I,0] = 0.0
158 shear.angvel[I,1] = 0.0
159 shear.angvel[I,2] = 0.0
160 shear.vel[I,0] = 0.0 # x-dim
161 shear.vel[I,1] = 0.0 # y-dim
162
163 # Copy bottom fixed particles to top
164 z_offset = z_max-z_min
165 shearvel = (z_max-z_min)*shear_strain_rate
166 for i in I[0]:
167 x = shear.x[i,:] + numpy.array([0.0, 0.0, z_offset])
168 vel = numpy.array([shearvel, 0.0, 0.0])
169 shear.addParticle(x = x, radius = shear.radius[i], fixvel = 1,
170 vel = vel)
171
172 # Set wall viscosities to zero
173 shear.gamma_wn[0] = 0.0
174 shear.gamma_wt[0] = 0.0
175
176 # Set wall friction coefficients to zero
177 shear.mu_ws[0] = 0.0
178 shear.mu_wd[0] = 0.0
179
180 # Readjust top wall
181 shear.adjustUpperWall()
182
183 # Set duration of simulation
184 #shear.initTemporal(total = 20.0)
185 shear.initTemporal(total = 400.0)
186
187 if (shearing == True):
188
189 # Run sphere
190 shear.run(dry=True)
191 shear.run()
192
193 if (plots == True):
194 # Make a graph of energies
195 shear.visualize('energy', savefig=True, outformat='png')
196 shear.visualize('shear', savefig=True, outformat='png')
197
198 if (rendering == True):
199 # Render images with raytracer
200 shear.render(method = 'pres', max_val = 2.0*devs, verbose = False)