world.py - sphere - GPU-based 3D discrete element method algorithm with optional fluid coupling
HTML git clone git://src.adamsgaard.dk/sphere
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
world.py (47043B)
---
1 import math
2 import numpy
3 from .common import FontProperties, V_sphere, plt, py_mpl
4
5
6 class SimSetup:
7 'Particle, grid, boundary and material setup for sim objects.'
8
9 def generateRadii(self, psd='logn', mean=440e-6, variance=8.8e-9,
10 histogram=False):
11 '''
12 Draw random particle radii from a selected probability distribution.
13 The larger the variance of radii is, the slower the computations will
14 run. The reason is two-fold: The smallest particle dictates the time
15 step length, where smaller particles cause shorter time steps. At the
16 same time, the largest particle determines the sorting cell size, where
17 larger particles cause larger cells. Larger cells are likely to contain
18 more particles, causing more contact checks.
19
20 :param psd: The particle side distribution. One possible value is
21 ``logn``, which is a log-normal probability distribution, suitable
22 for approximating well-sorted, coarse sediments. The other possible
23 value is ``uni``, which is a uniform distribution from
24 ``mean - variance`` to ``mean + variance``.
25 :type psd: str
26 :param mean: The mean radius [m] (default=440e-6 m)
27 :type mean: float
28 :param variance: The variance in the probability distribution
29 [m].
30 :type variance: float
31
32 See also: :func:`generateBimodalRadii()`.
33 '''
34
35 if psd == 'logn': # Log-normal probability distribution
36 mu = math.log((mean**2)/math.sqrt(variance+mean**2))
37 sigma = math.sqrt(math.log(variance/(mean**2)+1))
38 self.radius = numpy.random.lognormal(mu, sigma, self.np)
39 elif psd == 'uni': # Uniform distribution
40 radius_min = mean - variance
41 radius_max = mean + variance
42 self.radius = numpy.random.uniform(radius_min, radius_max, self.np)
43 else:
44 raise Exception('Particle size distribution type not understood ('
45 + str(psd) + '). '
46 + 'Valid values are \'uni\' or \'logn\'')
47
48 # Show radii as histogram
49 if histogram and py_mpl:
50 fig = plt.figure(figsize=(8, 8))
51 figtitle = 'Particle size distribution, {0} particles'\
52 .format(self.np)
53 fig.text(0.5, 0.95, figtitle, horizontalalignment='center',
54 fontproperties=FontProperties(size=18))
55 bins = 20
56
57 # Create histogram
58 plt.hist(self.radius, bins)
59
60 # Plot
61 plt.xlabel('Radii [m]')
62 plt.ylabel('Count')
63 plt.axis('tight')
64 fig.savefig(self.sid + '-psd.png')
65 fig.clf()
66
67 def generateBimodalRadii(self, r_small=0.005, r_large=0.05, ratio=0.2,
68 verbose=True):
69 '''
70 Draw random radii from two distinct sizes.
71
72 :param r_small: Radii of small population [m], in ]0;r_large[
73 :type r_small: float
74 :param r_large: Radii of large population [m], in ]r_small;inf[
75 :type r_large: float
76 :param ratio: Approximate volumetric ratio between the two
77 populations (large/small).
78 :type ratio: float
79
80 See also: :func:`generateRadii()`.
81 '''
82 if r_small >= r_large:
83 raise Exception("r_large should be larger than r_small")
84
85 V_small = V_sphere(r_small)
86 V_large = V_sphere(r_large)
87 nlarge = int(V_small/V_large * ratio * self.np) # ignore void volume
88
89 self.radius[:] = r_small
90 self.radius[0:nlarge] = r_large
91 numpy.random.shuffle(self.radius)
92
93 # Test volumetric ratio
94 V_small_total = V_small * (self.np - nlarge)
95 V_large_total = V_large * nlarge
96 if abs(V_large_total/V_small_total - ratio) > 1.0e5:
97 raise Exception("Volumetric ratio seems wrong")
98
99 if verbose:
100 print("generateBimodalRadii created " + str(nlarge)
101 + " large particles, and " + str(self.np - nlarge)
102 + " small")
103
104 def checkerboardColors(self, nx=6, ny=6, nz=6):
105 '''
106 Assign checkerboard color values to the particles in an orthogonal grid.
107
108 :param nx: Number of color values along the x axis
109 :type nx: int
110 :param ny: Number of color values along the y ayis
111 :type ny: int
112 :param nz: Number of color values along the z azis
113 :type nz: int
114 '''
115 x_min = numpy.min(self.x[:, 0])
116 x_max = numpy.max(self.x[:, 0])
117 y_min = numpy.min(self.x[:, 1])
118 y_max = numpy.max(self.x[:, 1])
119 z_min = numpy.min(self.x[:, 2])
120 z_max = numpy.max(self.x[:, 2])
121 for i in numpy.arange(self.np):
122 ix = numpy.floor((self.x[i, 0] - x_min)/(x_max/nx))
123 iy = numpy.floor((self.x[i, 1] - y_min)/(y_max/ny))
124 iz = numpy.floor((self.x[i, 2] - z_min)/(z_max/nz))
125 self.color[i] = (-1)**ix + (-1)**iy + (-1)**iz
126
127 def contactModel(self, contactmodel):
128 '''
129 Define which contact model to use for the tangential component of
130 particle-particle interactions. The elastic-viscous-frictional contact
131 model (2) is considered to be the most realistic contact model, while
132 the viscous-frictional contact model is significantly faster.
133
134 :param contactmodel: The type of tangential contact model to use
135 (visco-frictional=1, elasto-visco-frictional=2)
136 :type contactmodel: int
137 '''
138 self.contactmodel[0] = contactmodel
139
140 def wall0iz(self):
141 '''
142 Returns the cell index of wall 0 along z.
143
144 :returns: z cell index
145 :return type: int
146 '''
147 if self.nw > 0:
148 return int(self.w_x[0]/(self.L[2]/self.num[2]))
149 else:
150 raise Exception('No dynamic top wall present!')
151
152 def normalBoundariesXY(self):
153 '''
154 Set the x and y boundary conditions to be static walls.
155
156 See also :func:`periodicBoundariesXY()` and
157 :func:`periodicBoundariesX()`
158 '''
159 self.periodic[0] = 0
160
161 def periodicBoundariesXY(self):
162 '''
163 Set the x and y boundary conditions to be periodic.
164
165 See also :func:`normalBoundariesXY()` and
166 :func:`periodicBoundariesX()`
167 '''
168 self.periodic[0] = 1
169
170 def periodicBoundariesX(self):
171 '''
172 Set the x boundary conditions to be periodic.
173
174 See also :func:`normalBoundariesXY()` and
175 :func:`periodicBoundariesXY()`
176 '''
177 self.periodic[0] = 2
178
179 def adaptiveGrid(self):
180 '''
181 Set the height of the fluid grid to automatically readjust to the
182 height of the granular assemblage, as dictated by the position of the
183 top wall. This will readjust `self.L[2]` during the simulation to
184 equal the position of the top wall `self.w_x[0]`.
185
186 See also :func:`staticGrid()`
187 '''
188 self.adaptive[0] = 1
189
190 def staticGrid(self):
191 '''
192 Set the height of the fluid grid to be constant as set in `self.L[2]`.
193
194 See also :func:`adaptiveGrid()`
195 '''
196 self.adaptive[0] = 0
197
198 def initRandomPos(self, gridnum=numpy.array([12, 12, 36]), dx=-1.0):
199 '''
200 Initialize particle positions in completely random configuration. Radii
201 *must* be set beforehand. If the x and y boundaries are set as periodic,
202 the particle centers will be placed all the way to the edge. On regular,
203 non-periodic boundaries, the particles are restrained at the edges to
204 make space for their radii within the bounding box.
205
206 :param gridnum: The number of sorting cells in each spatial direction
207 (default=[12, 12, 36])
208 :type gridnum: numpy.array
209 :param dx: The cell width in any direction. If the default value is used
210 (-1), the cell width is calculated to fit the largest particle.
211 :type dx: float
212 '''
213
214 # Calculate cells in grid
215 self.num = gridnum
216 r_max = numpy.max(self.radius)
217
218 # Cell configuration
219 if dx > 0.0:
220 cellsize = dx
221 else:
222 cellsize = 2.1 * numpy.amax(self.radius)
223
224 # World size
225 self.L = self.num * cellsize
226
227 # Particle positions randomly distributed without overlap
228 for i in range(self.np):
229 overlaps = True
230 while overlaps:
231 overlaps = False
232
233 # Draw random position
234 for d in range(self.nd):
235 self.x[i, d] = (self.L[d] - self.origo[d] - 2*r_max) \
236 * numpy.random.random_sample() \
237 + self.origo[d] + r_max
238
239 # Check other particles for overlaps
240 for j in range(i-1):
241 delta = self.x[i] - self.x[j]
242 delta_len = math.sqrt(numpy.dot(delta, delta)) \
243 - (self.radius[i] + self.radius[j])
244 if delta_len < 0.0:
245 overlaps = True
246 print("\rFinding non-overlapping particle positions, "
247 + "{0} % complete".format(numpy.ceil(i/self.np*100)))
248
249 # Print newline
250 print()
251
252 def defineWorldBoundaries(self, L, origo=[0.0, 0.0, 0.0], dx=-1):
253 '''
254 Set the boundaries of the world. Particles will only be able to interact
255 within this domain. With dynamic walls, allow space for expansions.
256 *Important*: The particle radii have to be set beforehand. The world
257 edges act as static walls.
258
259 :param L: The upper boundary of the domain [m]
260 :type L: numpy.array
261 :param origo: The lower boundary of the domain [m]. Negative values
262 won't work. Default=[0.0, 0.0, 0.0].
263 :type origo: numpy.array
264 :param dx: The cell width in any direction. If the default value is used
265 (-1), the cell width is calculated to fit the largest particle.
266 :type dx: float
267 '''
268
269 # Cell configuration
270 if dx > 0.0:
271 cellsize_min = dx
272 else:
273 if self.np < 1:
274 raise Exception('Error: You need to define dx in ' +
275 'defineWorldBoundaries if there are no ' +
276 'particles in the simulation.')
277 cellsize_min = 2.1 * numpy.amax(self.radius)
278
279 # Lower boundary of the sorting grid
280 self.origo[:] = origo[:]
281
282 # Upper boundary of the sorting grid
283 self.L[:] = L[:]
284
285 # Adjust the number of sorting cells along each axis to fit the largest
286 # particle size and the world size
287 self.num[0] = numpy.ceil((self.L[0]-self.origo[0])/cellsize_min)
288 self.num[1] = numpy.ceil((self.L[1]-self.origo[1])/cellsize_min)
289 self.num[2] = numpy.ceil((self.L[2]-self.origo[2])/cellsize_min)
290
291 #if (self.num.any() < 4):
292 #if (self.num[0] < 4 or self.num[1] < 4 or self.num[2] < 4):
293 if self.num[0] < 3 or self.num[1] < 3 or self.num[2] < 3:
294 raise Exception("Error: The grid must be at least 3 cells in each "
295 + "direction\nGrid: x={}, y={}, z={}\n"
296 .format(self.num[0], self.num[1], self.num[2])
297 + "Please increase the world size.")
298
299 def initGrid(self, dx=-1):
300 '''
301 Initialize grid suitable for the particle positions set previously.
302 The margin parameter adjusts the distance (in no. of max. radii)
303 from the particle boundaries.
304 *Important*: The particle radii have to be set beforehand if the cell
305 width isn't specified by `dx`.
306
307 :param dx: The cell width in any direction. If the default value is used
308 (-1), the cell width is calculated to fit the largest particle.
309 :type dx: float
310 '''
311
312 # Cell configuration
313 if dx > 0.0:
314 cellsize_min = dx
315 else:
316 cellsize_min = 2.1 * numpy.amax(self.radius)
317 self.num[0] = numpy.ceil((self.L[0]-self.origo[0])/cellsize_min)
318 self.num[1] = numpy.ceil((self.L[1]-self.origo[1])/cellsize_min)
319 self.num[2] = numpy.ceil((self.L[2]-self.origo[2])/cellsize_min)
320
321 if self.num[0] < 4 or self.num[1] < 4 or self.num[2] < 4:
322 raise Exception("Error: The grid must be at least 3 cells in each "
323 + "direction\nGrid: x={}, y={}, z={}"
324 .format(self.num[0], self.num[1], self.num[2]))
325
326 # Put upper wall at top boundary
327 if self.nw > 0:
328 self.w_x[0] = self.L[0]
329
330 def initGridAndWorldsize(self, margin=2.0):
331 '''
332 Initialize grid suitable for the particle positions set previously.
333 The margin parameter adjusts the distance (in no. of max. radii)
334 from the particle boundaries. If the upper wall is dynamic, it is placed
335 at the top boundary of the world.
336
337 :param margin: Distance to world boundary in no. of max. particle radii
338 :type margin: float
339 '''
340
341 # Cell configuration
342 r_max = numpy.amax(self.radius)
343
344 # Max. and min. coordinates of world
345 self.origo = numpy.array([numpy.amin(self.x[:, 0] - self.radius[:]),
346 numpy.amin(self.x[:, 1] - self.radius[:]),
347 numpy.amin(self.x[:, 2] - self.radius[:])]) \
348 - margin*r_max
349 self.L = numpy.array([numpy.amax(self.x[:, 0] + self.radius[:]),
350 numpy.amax(self.x[:, 1] + self.radius[:]),
351 numpy.amax(self.x[:, 2] + self.radius[:])]) \
352 + margin*r_max
353
354 cellsize_min = 2.1 * r_max
355 self.num[0] = numpy.ceil((self.L[0]-self.origo[0])/cellsize_min)
356 self.num[1] = numpy.ceil((self.L[1]-self.origo[1])/cellsize_min)
357 self.num[2] = numpy.ceil((self.L[2]-self.origo[2])/cellsize_min)
358
359 if self.num[0] < 4 or self.num[1] < 4 or self.num[2] < 4:
360 raise Exception("Error: The grid must be at least 3 cells in each "
361 + "direction, num=" + str(self.num))
362
363 # Put upper wall at top boundary
364 if self.nw > 0:
365 self.w_x[0] = self.L[0]
366
367 def initGridPos(self, gridnum=numpy.array([12, 12, 36])):
368 '''
369 Initialize particle positions in loose, cubic configuration.
370 ``gridnum`` is the number of cells in the x, y and z directions.
371 *Important*: The particle radii and the boundary conditions (periodic or
372 not) for the x and y boundaries have to be set beforehand.
373
374 :param gridnum: The number of particles in x, y and z directions
375 :type gridnum: numpy.array
376 '''
377
378 # Calculate cells in grid
379 self.num = numpy.asarray(gridnum)
380
381 # World size
382 r_max = numpy.amax(self.radius)
383 cellsize = 2.1 * r_max
384 self.L = self.num * cellsize
385
386 # Check whether there are enough grid cells
387 if (self.num[0]*self.num[1]*self.num[2]-(2**3)) < self.np:
388 print("Error! The grid is not sufficiently large.")
389 raise NameError('Error! The grid is not sufficiently large.')
390
391 gridpos = numpy.zeros(self.nd, dtype=numpy.uint32)
392
393 # Make sure grid is sufficiently large if every second level is moved
394 if self.periodic[0] == 1:
395 self.num[0] -= 1
396 self.num[1] -= 1
397
398 # Check whether there are enough grid cells
399 if (self.num[0]*self.num[1]*self.num[2]-(2*3*3)) < self.np:
400 print("Error! The grid is not sufficiently large.")
401 raise NameError('Error! The grid is not sufficiently large.')
402
403 # Particle positions randomly distributed without overlap
404 for i in range(self.np):
405
406 # Find position in 3d mesh from linear index
407 gridpos[0] = (i % (self.num[0]))
408 gridpos[1] = numpy.floor(i/(self.num[0])) % (self.num[0])
409 gridpos[2] = numpy.floor(i/((self.num[0])*(self.num[1]))) #\
410 #% ((self.num[0])*(self.num[1]))
411
412 for d in range(self.nd):
413 self.x[i, d] = gridpos[d] * cellsize + 0.5*cellsize
414
415 # Allow pushing every 2.nd level out of lateral boundaries
416 if self.periodic[0] == 1:
417 # Offset every second level
418 if gridpos[2] % 2:
419 self.x[i, 0] += 0.5*cellsize
420 self.x[i, 1] += 0.5*cellsize
421
422 # Readjust grid to correct size
423 if self.periodic[0] == 1:
424 self.num[0] += 1
425 self.num[1] += 1
426
427 def initRandomGridPos(self, gridnum=numpy.array([12, 12, 32]),
428 padding=2.1):
429 '''
430 Initialize particle positions in loose, cubic configuration with some
431 variance. ``gridnum`` is the number of cells in the x, y and z
432 directions. *Important*: The particle radii and the boundary conditions
433 (periodic or not) for the x and y boundaries have to be set beforehand.
434 The world size and grid height (in the z direction) is readjusted to fit
435 the particle positions.
436
437 :param gridnum: The number of particles in x, y and z directions
438 :type gridnum: numpy.array
439 :param padding: Increase distance between particles in x, y and z
440 directions with this multiplier. Large values create more random
441 packings.
442 :type padding: float
443 '''
444
445 # Calculate cells in grid
446 coarsegrid = numpy.floor(numpy.asarray(gridnum)/2)
447
448 # World size
449 r_max = numpy.amax(self.radius)
450
451 # Cells in grid 2*size to make space for random offset
452 cellsize = padding * r_max * 2
453
454 # Check whether there are enough grid cells
455 if ((coarsegrid[0]-1)*(coarsegrid[1]-1)*(coarsegrid[2]-1)) < self.np:
456 print("Error! The grid is not sufficiently large.")
457 raise NameError('Error! The grid is not sufficiently large.')
458
459 gridpos = numpy.zeros(self.nd, dtype=numpy.uint32)
460
461 # Particle positions randomly distributed without overlap
462 for i in range(self.np):
463
464 # Find position in 3d mesh from linear index
465 gridpos[0] = (i % (coarsegrid[0]))
466 gridpos[1] = numpy.floor(i/(coarsegrid[0]))%(coarsegrid[1]) # Thanks Horacio!
467 gridpos[2] = numpy.floor(i/((coarsegrid[0])*(coarsegrid[1])))
468
469 # Place particles in grid structure, and randomly adjust the
470 # positions within the oversized cells (uniform distribution)
471 for d in range(self.nd):
472 r = self.radius[i]*1.05
473 self.x[i, d] = gridpos[d] * cellsize \
474 + ((cellsize-r) - r) \
475 * numpy.random.random_sample() + r
476
477 # Calculate new grid with cell size equal to max. particle diameter
478 x_max = numpy.max(self.x[:, 0] + self.radius)
479 y_max = numpy.max(self.x[:, 1] + self.radius)
480 z_max = numpy.max(self.x[:, 2] + self.radius)
481
482 # Adjust size of world
483 self.num[0] = numpy.ceil(x_max/cellsize)
484 self.num[1] = numpy.ceil(y_max/cellsize)
485 self.num[2] = numpy.ceil(z_max/cellsize)
486 self.L = self.num * cellsize
487
488 def createBondPair(self, i, j, spacing=-0.1):
489 '''
490 Bond particles i and j. Particle j is moved adjacent to particle i,
491 and oriented randomly.
492
493 :param i: Index of first particle in bond
494 :type i: int
495 :param j: Index of second particle in bond
496 :type j: int
497 :param spacing: The inter-particle distance prescribed. Positive
498 values result in a inter-particle distance, negative equal an
499 overlap. The value is relative to the sum of the two radii.
500 :type spacing: float
501 '''
502
503 x_i = self.x[i]
504 r_i = self.radius[i]
505 r_j = self.radius[j]
506 dist_ij = (r_i + r_j)*(1.0 + spacing)
507
508 dazi = numpy.random.rand(1) * 360.0 # azimuth
509 azi = numpy.radians(dazi)
510 dang = numpy.random.rand(1) * 180.0 - 90.0 # angle
511 ang = numpy.radians(dang)
512
513 x_j = numpy.copy(x_i)
514 x_j[0] = x_j[0] + dist_ij * numpy.cos(azi) * numpy.cos(ang)
515 x_j[1] = x_j[1] + dist_ij * numpy.sin(azi) * numpy.cos(ang)
516 x_j[2] = x_j[2] + dist_ij * numpy.sin(ang) * numpy.cos(azi)
517 self.x[j] = x_j
518
519 if self.x[j, 0] < self.origo[0]:
520 self.x[j, 0] += x_i[0] - x_j[0]
521 if self.x[j, 1] < self.origo[1]:
522 self.x[j, 1] += x_i[1] - x_j[1]
523 if self.x[j, 2] < self.origo[2]:
524 self.x[j, 2] += x_i[2] - x_j[2]
525
526 if self.x[j, 0] > self.L[0]:
527 self.x[j, 0] -= abs(x_j[0] - x_i[0])
528 if self.x[j, 1] > self.L[1]:
529 self.x[j, 1] -= abs(x_j[1] - x_i[1])
530 if self.x[j, 2] > self.L[2]:
531 self.x[j, 2] -= abs(x_j[2] - x_i[2])
532
533 self.bond(i, j) # register bond
534
535 # Check that the spacing is correct
536 x_ij = self.x[i] - self.x[j]
537 x_ij_length = numpy.sqrt(x_ij.dot(x_ij))
538 if (x_ij_length - dist_ij) > dist_ij*0.01:
539 print(x_i); print(r_i)
540 print(x_j); print(r_j)
541 print(x_ij_length); print(dist_ij)
542 raise Exception("Error, something went wrong in createBondPair")
543
544 def randomBondPairs(self, ratio=0.3, spacing=-0.1):
545 '''
546 Bond an amount of particles in two-particle clusters. The particles
547 should be initialized beforehand. Note: The actual number of bonds is
548 likely to be somewhat smaller than specified, due to the random
549 selection algorithm.
550
551 :param ratio: The amount of particles to bond, values in ]0.0;1.0]
552 :type ratio: float
553 :param spacing: The distance relative to the sum of radii between bonded
554 particles, neg. values denote an overlap. Values in ]0.0,inf[.
555 :type spacing: float
556 '''
557
558 bondparticles = numpy.unique(numpy.random.random_integers(0, high=self.np-1,
559 size=int(self.np*ratio)))
560 if bondparticles.size % 2 > 0:
561 bondparticles = bondparticles[:-1].copy()
562 bondparticles = bondparticles.reshape(int(bondparticles.size/2),
563 2).copy()
564
565 for n in numpy.arange(bondparticles.shape[0]):
566 self.createBondPair(bondparticles[n, 0], bondparticles[n, 1],
567 spacing)
568
569 def zeroKinematics(self):
570 '''
571 Zero all kinematic parameters of the particles. This function is useful
572 when output from one simulation is reused in another simulation.
573 '''
574
575 self.force = numpy.zeros((self.np, self.nd))
576 self.torque = numpy.zeros((self.np, self.nd))
577 self.vel = numpy.zeros(self.np*self.nd, dtype=numpy.float64)\
578 .reshape(self.np, self.nd)
579 self.angvel = numpy.zeros(self.np*self.nd, dtype=numpy.float64)\
580 .reshape(self.np, self.nd)
581 self.angpos = numpy.zeros(self.np*self.nd, dtype=numpy.float64)\
582 .reshape(self.np, self.nd)
583 self.es = numpy.zeros(self.np, dtype=numpy.float64)
584 self.ev = numpy.zeros(self.np, dtype=numpy.float64)
585 self.xyzsum = numpy.zeros(self.np*3, dtype=numpy.float64).reshape(self.np, 3)
586
587 def adjustUpperWall(self, z_adjust=1.1):
588 '''
589 Included for legacy purposes, calls :func:`adjustWall()` with ``idx=0``.
590
591 :param z_adjust: Increase the world and grid size by this amount to
592 allow for wall movement.
593 :type z_adjust: float
594 '''
595
596 # Initialize upper wall
597 self.nw = 1
598 self.wmode = numpy.zeros(1) # fixed BC
599 self.w_n = numpy.zeros(self.nw*self.nd, dtype=numpy.float64)\
600 .reshape(self.nw, self.nd)
601 self.w_n[0, 2] = -1.0
602 self.w_vel = numpy.zeros(1)
603 self.w_force = numpy.zeros(1)
604 self.w_sigma0 = numpy.zeros(1)
605
606 self.w_x = numpy.zeros(1)
607 self.w_m = numpy.zeros(1)
608 self.adjustWall(idx=0, adjust=z_adjust)
609
610 def adjustWall(self, idx, adjust=1.1):
611 '''
612 Adjust grid and dynamic wall to max. particle position. The wall
613 thickness will by standard equal the maximum particle diameter. The
614 density equals the particle density, and the wall size is equal to the
615 width and depth of the simulation domain (`self.L[0]` and `self.L[1]`).
616
617 :param: idx: The wall to adjust. 0=+z, upper wall (default), 1=-x,
618 left wall, 2=+x, right wall, 3=-y, front wall, 4=+y, back
619 wall.
620 :type idx: int
621 :param adjust: Increase the world and grid size by this amount to
622 allow for wall movement.
623 :type adjust: float
624 '''
625
626 if idx == 0:
627 dim = 2
628 elif idx == 1 or idx == 2:
629 dim = 0
630 elif idx == 3 or idx == 4:
631 dim = 1
632 else:
633 print("adjustWall: idx value not understood")
634
635 xmin = numpy.min(self.x[:, dim] - self.radius)
636 xmax = numpy.max(self.x[:, dim] + self.radius)
637
638 cellsize = self.L[0] / self.num[0]
639 self.num[dim] = numpy.ceil(((xmax-xmin)*adjust + xmin)/cellsize)
640 self.L[dim] = (xmax-xmin)*adjust + xmin
641
642 # Initialize upper wall
643 if idx == 0 or idx == 1 or idx == 3:
644 self.w_x[idx] = xmax
645 else:
646 self.w_x[idx] = xmin
647 self.w_m[idx] = self.totalMass()
648
649 def consolidate(self, normal_stress=10e3):
650 '''
651 Setup consolidation experiment. Specify the upper wall normal stress in
652 Pascal, default value is 10 kPa.
653
654 :param normal_stress: The normal stress to apply from the upper wall
655 :type normal_stress: float
656 '''
657
658 self.nw = 1
659
660 if normal_stress <= 0.0:
661 raise Exception('consolidate() error: The normal stress should be '
662 'a positive value, but is ' + str(normal_stress) +
663 ' Pa')
664
665 # Zero the kinematics of all particles
666 self.zeroKinematics()
667
668 # Adjust grid and placement of upper wall
669 self.adjustUpperWall()
670
671 # Set the top wall BC to a value of normal stress
672 self.wmode = numpy.array([1])
673 self.w_sigma0 = numpy.ones(1) * normal_stress
674
675 # Set top wall to a certain mass corresponding to the selected normal
676 # stress
677 #self.w_sigma0 = numpy.zeros(1)
678 #self.w_m[0] = numpy.abs(normal_stress*self.L[0]*self.L[1]/self.g[2])
679 self.w_m[0] = self.totalMass()
680
681 def uniaxialStrainRate(self, wvel=-0.001):
682 '''
683 Setup consolidation experiment. Specify the upper wall velocity in m/s,
684 default value is -0.001 m/s (i.e. downwards).
685
686 :param wvel: Upper wall velocity. Negative values mean that the wall
687 moves downwards.
688 :type wvel: float
689 '''
690
691 # zero kinematics
692 self.zeroKinematics()
693
694 # Initialize upper wall
695 self.adjustUpperWall()
696 self.wmode = numpy.array([2]) # strain rate BC
697 self.w_vel = numpy.array([wvel])
698
699 def triaxial(self, wvel=-0.001, normal_stress=10.0e3):
700 '''
701 Setup triaxial experiment. The upper wall is moved at a fixed velocity
702 in m/s, default values is -0.001 m/s (i.e. downwards). The side walls
703 are exerting a defined normal stress.
704
705 :param wvel: Upper wall velocity. Negative values mean that the wall
706 moves downwards.
707 :type wvel: float
708 :param normal_stress: The normal stress to apply from the upper wall.
709 :type normal_stress: float
710 '''
711
712 # zero kinematics
713 self.zeroKinematics()
714
715 # Initialize walls
716 self.nw = 5 # five dynamic walls
717 self.wmode = numpy.array([2, 1, 1, 1, 1]) # BCs (vel, stress, stress, ...)
718 self.w_vel = numpy.array([1, 0, 0, 0, 0]) * wvel
719 self.w_sigma0 = numpy.array([0, 1, 1, 1, 1]) * normal_stress
720 self.w_n = numpy.array(([0, 0, -1], [-1, 0, 0],
721 [1, 0, 0], [0, -1, 0], [0, 1, 0]),
722 dtype=numpy.float64)
723 self.w_x = numpy.zeros(5)
724 self.w_m = numpy.zeros(5)
725 self.w_force = numpy.zeros(5)
726 for i in range(5):
727 self.adjustWall(idx=i)
728
729 def shear(self, shear_strain_rate=1.0, shear_stress=False):
730 '''
731 Setup shear experiment either by a constant shear rate or a constant
732 shear stress. The shear strain rate is the shear velocity divided by
733 the initial height per second. The shear movement is along the positive
734 x axis. The function zeroes the tangential wall viscosity (gamma_wt) and
735 the wall friction coefficients (mu_ws, mu_wn).
736
737 :param shear_strain_rate: The shear strain rate [-] to use if
738 shear_stress isn't False.
739 :type shear_strain_rate: float
740 :param shear_stress: The shear stress value to use [Pa].
741 :type shear_stress: float or bool
742 '''
743
744 self.nw = 1
745
746 # Find lowest and heighest point
747 z_min = numpy.min(self.x[:, 2] - self.radius)
748 z_max = numpy.max(self.x[:, 2] + self.radius)
749
750 # the grid cell size is equal to the max. particle diameter
751 cellsize = self.L[0] / self.num[0]
752
753 # make grid one cell heigher to allow dilation
754 self.num[2] += 1
755 self.L[2] = self.num[2] * cellsize
756
757 # zero kinematics
758 self.zeroKinematics()
759
760 # Adjust grid and placement of upper wall
761 self.wmode = numpy.array([1])
762
763 # Fix horizontal velocity to 0.0 of lowermost particles
764 d_max_below = numpy.max(self.radius[numpy.nonzero(self.x[:, 2] <
765 (z_max-z_min)*0.3)])*2.0
766 I = numpy.nonzero(self.x[:, 2] < (z_min + d_max_below))
767 self.fixvel[I] = 1
768 self.angvel[I, 0] = 0.0
769 self.angvel[I, 1] = 0.0
770 self.angvel[I, 2] = 0.0
771 self.vel[I, 0] = 0.0 # x-dim
772 self.vel[I, 1] = 0.0 # y-dim
773 self.color[I] = -1
774
775 # Fix horizontal velocity to specific value of uppermost particles
776 d_max_top = numpy.max(self.radius[numpy.nonzero(self.x[:, 2] >
777 (z_max-z_min)*0.7)])*2.0
778 I = numpy.nonzero(self.x[:, 2] > (z_max - d_max_top))
779 self.fixvel[I] = 1
780 self.angvel[I, 0] = 0.0
781 self.angvel[I, 1] = 0.0
782 self.angvel[I, 2] = 0.0
783 if not shear_stress:
784 self.vel[I, 0] = (z_max-z_min)*shear_strain_rate
785 else:
786 self.vel[I, 0] = 0.0
787 self.wmode[0] = 3
788 self.w_tau_x[0] = float(shear_stress)
789 self.vel[I, 1] = 0.0 # y-dim
790 self.color[I] = -1
791
792 # Set wall tangential viscosity to zero
793 self.gamma_wt[0] = 0.0
794
795 # Set wall friction coefficients to zero
796 self.mu_ws[0] = 0.0
797 self.mu_wd[0] = 0.0
798
799 def initTemporal(self, total, current=0.0, file_dt=0.05, step_count=0,
800 dt=-1, epsilon=0.01):
801 '''
802 Set temporal parameters for the simulation. *Important*: Particle radii,
803 physical parameters, and the optional fluid grid need to be set prior to
804 these if the computational time step (dt) isn't set explicitly. If the
805 parameter `dt` is the default value (-1), the function will estimate the
806 best time step length. The value of the computational time step for the
807 DEM is checked for stability in the CFD solution if fluid simulation is
808 included.
809
810 :param total: The time at which to end the simulation [s]
811 :type total: float
812 :param current: The current time [s] (default=0.0 s)
813 :type total: float
814 :param file_dt: The interval between output files [s] (default=0.05 s)
815 :type total: float
816 :step_count: The number of the first output file (default=0)
817 :type step_count: int
818 :param dt: The computational time step length [s]
819 :type total: float
820 :param epsilon: Time step multiplier (default=0.01)
821 :type epsilon: float
822 '''
823
824 if dt > 0.0:
825 self.time_dt[0] = dt
826 if self.np > 0:
827 print("Warning: Manually specifying the time step length when "
828 "simulating particles may produce instabilities.")
829
830 elif self.np > 0:
831
832 r_min = numpy.min(self.radius)
833 m_min = self.rho[0] * 4.0/3.0*numpy.pi*r_min**3
834
835 if self.E > 0.001:
836 k_max = numpy.max(numpy.pi/2.0*self.E*self.radius)
837 else:
838 k_max = numpy.max([self.k_n[:], self.k_t[:]])
839
840 # Radjaii et al 2011
841 self.time_dt[0] = epsilon/(numpy.sqrt(k_max/m_min))
842
843 # Zhang and Campbell, 1992
844 #self.time_dt[0] = 0.075*math.sqrt(m_min/k_max)
845
846 # Computational time step (O'Sullivan et al, 2003)
847 #self.time_dt[0] = 0.17*math.sqrt(m_min/k_max)
848
849 elif not self.fluid:
850 raise Exception('Error: Could not automatically set a time step.')
851
852 # Check numerical stability of the fluid phase, by criteria derived
853 # by von Neumann stability analysis of the diffusion and advection
854 # terms
855 if self.fluid:
856 fluid_time_dt = self.largestFluidTimeStep()
857 self.time_dt[0] = numpy.min([fluid_time_dt, self.time_dt[0]])
858
859 # Time at start
860 self.time_current[0] = current
861 self.time_total[0] = total
862 self.time_file_dt[0] = file_dt
863 self.time_step_count[0] = step_count
864
865 def currentTime(self, value=-1):
866 '''
867 Get or set the current time. If called without arguments the current
868 time is returned. If a new time is passed in the 'value' argument, the
869 time is written to the object.
870
871 :param value: The new current time
872 :type value: float
873
874 :returns: The current time
875 :return type: float
876 '''
877 if value != -1:
878 self.time_current[0] = value
879 else:
880 return self.time_current[0]
881
882 def defaultParams(self, mu_s=0.5, mu_d=0.5, mu_r=0.0, rho=2600, k_n=1.16e9,
883 k_t=1.16e9, k_r=0, gamma_n=0.0, gamma_t=0.0, gamma_r=0.0,
884 gamma_wn=0.0, gamma_wt=0.0, capillaryCohesion=0):
885 '''
886 Initialize particle parameters to default values.
887
888 :param mu_s: The coefficient of static friction between particles [-]
889 :type mu_s: float
890 :param mu_d: The coefficient of dynamic friction between particles [-]
891 :type mu_d: float
892 :param rho: The density of the particle material [kg/(m^3)]
893 :type rho: float
894 :param k_n: The normal stiffness of the particles [N/m]
895 :type k_n: float
896 :param k_t: The tangential stiffness of the particles [N/m]
897 :type k_t: float
898 :param k_r: The rolling stiffness of the particles [N/rad] *Parameter
899 not used*
900 :type k_r: float
901 :param gamma_n: Particle-particle contact normal viscosity [Ns/m]
902 :type gamma_n: float
903 :param gamma_t: Particle-particle contact tangential viscosity [Ns/m]
904 :type gamma_t: float
905 :param gamma_r: Particle-particle contact rolling viscosity *Parameter
906 not used*
907 :type gamma_r: float
908 :param gamma_wn: Wall-particle contact normal viscosity [Ns/m]
909 :type gamma_wn: float
910 :param gamma_wt: Wall-particle contact tangential viscosity [Ns/m]
911 :type gamma_wt: float
912 :param capillaryCohesion: Enable particle-particle capillary cohesion
913 interaction model (0=no (default), 1=yes)
914 :type capillaryCohesion: int
915 '''
916
917 # Particle material density, kg/m^3
918 self.rho = numpy.ones(1, dtype=numpy.float64) * rho
919
920
921 ### Dry granular material parameters
922
923 # Contact normal elastic stiffness, N/m
924 self.k_n = numpy.ones(1, dtype=numpy.float64) * k_n
925
926 # Contact shear elastic stiffness (for contactmodel=2), N/m
927 self.k_t = numpy.ones(1, dtype=numpy.float64) * k_t
928
929 # Contact rolling elastic stiffness (for contactmodel=2), N/m
930 self.k_r = numpy.ones(1, dtype=numpy.float64) * k_r
931
932 # Contact normal viscosity. Critical damping: 2*sqrt(m*k_n).
933 # Normal force component elastic if nu=0.0.
934 #self.gamma_n=numpy.ones(self.np, dtype=numpy.float64) \
935 # * nu_frac * 2.0 * math.sqrt(4.0/3.0 * math.pi \
936 # * numpy.amin(self.radius)**3 \
937 # * self.rho[0] * self.k_n[0])
938 self.gamma_n = numpy.ones(1, dtype=numpy.float64) * gamma_n
939
940 # Contact shear viscosity, Ns/m
941 self.gamma_t = numpy.ones(1, dtype=numpy.float64) * gamma_t
942
943 # Contact rolling viscosity, Ns/m?
944 self.gamma_r = numpy.ones(1, dtype=numpy.float64) * gamma_r
945
946 # Contact static shear friction coefficient
947 #self.mu_s = numpy.ones(1, dtype=numpy.float64) * \
948 #numpy.tan(numpy.radians(ang_s))
949 self.mu_s = numpy.ones(1, dtype=numpy.float64) * mu_s
950
951 # Contact dynamic shear friction coefficient
952 #self.mu_d = numpy.ones(1, dtype=numpy.float64) * \
953 #numpy.tan(numpy.radians(ang_d))
954 self.mu_d = numpy.ones(1, dtype=numpy.float64) * mu_d
955
956 # Contact rolling friction coefficient
957 #self.mu_r = numpy.ones(1, dtype=numpy.float64) * \
958 #numpy.tan(numpy.radians(ang_r))
959 self.mu_r = numpy.ones(1, dtype=numpy.float64) * mu_r
960
961 # Wall viscosities
962 self.gamma_wn[0] = gamma_wn # normal
963 self.gamma_wt[0] = gamma_wt # sliding
964
965 # Wall friction coefficients
966 self.mu_ws = self.mu_s # static
967 self.mu_wd = self.mu_d # dynamic
968
969 ### Parameters related to capillary bonds
970
971 # Wettability, 0=perfect
972 theta = 0.0
973 if capillaryCohesion == 1:
974 # Prefactor
975 self.kappa[0] = 2.0 * math.pi * gamma_t * numpy.cos(theta)
976 self.V_b[0] = 1e-12 # Liquid volume at bond
977 else:
978 self.kappa[0] = 0.0 # Zero capillary force
979 self.V_b[0] = 0.0 # Zero liquid volume at bond
980
981 # Debonding distance
982 self.db[0] = (1.0 + theta/2.0) * self.V_b[0]**(1.0/3.0)
983
984 def setStiffnessNormal(self, k_n):
985 '''
986 Set the elastic stiffness (`k_n`) in the normal direction of the
987 contact.
988
989 :param k_n: The elastic stiffness coefficient [N/m]
990 :type k_n: float
991 '''
992 self.k_n[0] = k_n
993
994 def setStiffnessTangential(self, k_t):
995 '''
996 Set the elastic stiffness (`k_t`) in the tangential direction of the
997 contact.
998
999 :param k_t: The elastic stiffness coefficient [N/m]
1000 :type k_t: float
1001 '''
1002 self.k_t[0] = k_t
1003
1004 def setYoungsModulus(self, E):
1005 '''
1006 Set the elastic Young's modulus (`E`) for the contact model. This
1007 parameter is used over normal stiffness (`k_n`) and tangential
1008 stiffness (`k_t`) when its value is greater than zero. Using this
1009 parameter produces size-invariant behavior.
1010
1011 Example values are ~70e9 Pa for quartz,
1012 http://www.engineeringtoolbox.com/young-modulus-d_417.html
1013
1014 :param E: The elastic modulus [Pa]
1015 :type E: float
1016 '''
1017 self.E[0] = E
1018
1019 def setDampingNormal(self, gamma, over_damping=False):
1020 '''
1021 Set the dampening coefficient (gamma) in the normal direction of the
1022 particle-particle contact model. The function will print the fraction
1023 between the chosen damping and the critical damping value.
1024
1025 :param gamma: The viscous damping constant [N/(m/s)]
1026 :type gamma: float
1027 :param over_damping: Accept overdampening
1028 :type over_damping: boolean
1029
1030 See also: :func:`setDampingTangential(gamma)`
1031 '''
1032 self.gamma_n[0] = gamma
1033 critical_gamma = 2.0*numpy.sqrt(self.smallestMass()*self.k_n[0])
1034 damping_ratio = gamma/critical_gamma
1035 if damping_ratio < 1.0:
1036 print('Info: The system is under-dampened (ratio='
1037 + str(damping_ratio)
1038 + ') in the normal component. \nCritical damping='
1039 + str(critical_gamma) + '. This is ok.')
1040 elif damping_ratio > 1.0:
1041 if over_damping:
1042 print('Warning: The system is over-dampened (ratio='
1043 + str(damping_ratio) + ') in the normal component. '
1044 '\nCritical damping=' + str(critical_gamma) + '.')
1045 else:
1046 raise Exception('Warning: The system is over-dampened (ratio='
1047 + str(damping_ratio) + ') in the normal '
1048 'component.\n'
1049 'Call this function once more with '
1050 '`over_damping=True` if this is what you want.'
1051 '\nCritical damping=' + str(critical_gamma) +
1052 '.')
1053 else:
1054 print('Warning: The system is critically dampened (ratio=' +
1055 str(damping_ratio) + ') in the normal component. '
1056 '\nCritical damping=' + str(critical_gamma) + '.')
1057
1058 def setDampingTangential(self, gamma, over_damping=False):
1059 '''
1060 Set the dampening coefficient (gamma) in the tangential direction of the
1061 particle-particle contact model. The function will print the fraction
1062 between the chosen damping and the critical damping value.
1063
1064 :param gamma: The viscous damping constant [N/(m/s)]
1065 :type gamma: float
1066 :param over_damping: Accept overdampening
1067 :type over_damping: boolean
1068
1069 See also: :func:`setDampingNormal(gamma)`
1070 '''
1071 self.gamma_t[0] = gamma
1072 damping_ratio = gamma/(2.0*numpy.sqrt(self.smallestMass()*self.k_t[0]))
1073 if damping_ratio < 1.0:
1074 print('Info: The system is under-dampened (ratio='
1075 + str(damping_ratio)
1076 + ') in the tangential component. This is ok.')
1077 elif damping_ratio > 1.0:
1078 if over_damping:
1079 print('Warning: The system is over-dampened (ratio='
1080 + str(damping_ratio) + ') in the tangential component.')
1081 else:
1082 raise Exception('Warning: The system is over-dampened (ratio='
1083 + str(damping_ratio) + ') in the tangential '
1084 'component.\n'
1085 'Call this function once more with '
1086 '`over_damping=True` if this is what you want.')
1087 else:
1088 print('Warning: The system is critically dampened (ratio='
1089 + str(damping_ratio) + ') in the tangential component.')
1090
1091 def setStaticFriction(self, mu_s):
1092 '''
1093 Set the static friction coefficient for particle-particle interactions
1094 (`self.mu_s`). This value describes the resistance to a shearing motion
1095 while it is not happenind (contact tangential velocity zero).
1096
1097 :param mu_s: Value of the static friction coefficient, in [0;inf[.
1098 Usually between 0 and 1.
1099 :type mu_s: float
1100
1101 See also: :func:`setDynamicFriction(mu_d)`
1102 '''
1103 self.mu_s[0] = mu_s
1104
1105 def setDynamicFriction(self, mu_d):
1106 '''
1107 Set the dynamic friction coefficient for particle-particle interactions
1108 (`self.mu_d`). This value describes the resistance to a shearing motion
1109 while it is happening (contact tangential velocity larger than 0).
1110 Strain softening can be introduced by having a smaller dynamic
1111 frictional coefficient than the static fricion coefficient. Usually this
1112 value is identical to the static friction coefficient.
1113
1114 :param mu_d: Value of the dynamic friction coefficient, in [0;inf[.
1115 Usually between 0 and 1.
1116 :type mu_d: float
1117
1118 See also: :func:`setStaticFriction(mu_s)`
1119 '''
1120 self.mu_d[0] = mu_d
1121
1122 def scaleSize(self, factor):
1123 '''
1124 Scale the positions, linear velocities, forces, torques and radii of all
1125 particles and mobile walls.
1126
1127 :param factor: Spatial scaling factor ]0;inf[
1128 :type factor: float
1129 '''
1130 self.L *= factor
1131 self.x *= factor
1132 self.radius *= factor
1133 self.xyzsum *= factor
1134 self.vel *= factor
1135 self.force *= factor
1136 self.torque *= factor
1137 self.w_x *= factor
1138 self.w_m *= factor
1139 self.w_vel *= factor
1140 self.w_force *= factor
1141
1142 def bond(self, i, j):
1143 '''
1144 Create a bond between particles with index i and j
1145
1146 :param i: Index of first particle in bond
1147 :type i: int
1148 :param j: Index of second particle in bond
1149 :type j: int
1150 '''
1151
1152 self.lambda_bar[0] = 1.0 # Radius multiplier to parallel-bond radii
1153
1154 if not hasattr(self, 'bonds'):
1155 self.bonds = numpy.array([[i, j]], dtype=numpy.uint32)
1156 else:
1157 self.bonds = numpy.vstack((self.bonds, [i, j]))
1158
1159 if not hasattr(self, 'bonds_delta_n'):
1160 self.bonds_delta_n = numpy.array([0.0], dtype=numpy.uint32)
1161 else:
1162 #self.bonds_delta_n = numpy.vstack((self.bonds_delta_n, [0.0]))
1163 self.bonds_delta_n = numpy.append(self.bonds_delta_n, [0.0])
1164
1165 if not hasattr(self, 'bonds_delta_t'):
1166 self.bonds_delta_t = numpy.array([[0.0, 0.0, 0.0]], dtype=numpy.uint32)
1167 else:
1168 self.bonds_delta_t = numpy.vstack((self.bonds_delta_t,
1169 [0.0, 0.0, 0.0]))
1170
1171 if not hasattr(self, 'bonds_omega_n'):
1172 self.bonds_omega_n = numpy.array([0.0], dtype=numpy.uint32)
1173 else:
1174 #self.bonds_omega_n = numpy.vstack((self.bonds_omega_n, [0.0]))
1175 self.bonds_omega_n = numpy.append(self.bonds_omega_n, [0.0])
1176
1177 if not hasattr(self, 'bonds_omega_t'):
1178 self.bonds_omega_t = numpy.array([[0.0, 0.0, 0.0]],
1179 dtype=numpy.uint32)
1180 else:
1181 self.bonds_omega_t = numpy.vstack((self.bonds_omega_t,
1182 [0.0, 0.0, 0.0]))
1183
1184 # Increment the number of bonds with one
1185 self.nb0 += 1