tcapillary-cohesion2.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
---
tcapillary-cohesion2.py (3651B)
---
1 #!/usr/bin/env python
2
3 # This script simulates the effect of capillary cohesion on a sand pile put on a
4 # desk.
5
6 # start with
7 # $ python capillary-cohesion.py <DEVICE> <COHESION>
8 # where DEVICE specifies the index of the GPU (0 is the most common value).
9 # COHESION should have the value of 0 or 1. 0 denotes a dry simulation without
10 # cohesion, 1 denotes a wet simulation with capillary cohesion.
11 # GRAVITY toggles gravitational acceleration. Without it, the particles are
12 # placed in the middle of a volume. With it enabled, the particles are put on
13 # top of a flat wall.
14
15 import sphere
16 import numpy
17 import sys
18
19 device = sys.argv[1]
20 cohesion = sys.argv[2]
21
22 cube = sphere.sim('cube-init')
23 cube.readlast()
24 cube.adjustUpperWall(z_adjust=1.0)
25
26 # shrink particles to new mean radius and resize domain
27 r_mean = 0.001
28 r_mean_old = numpy.mean(cube.radius)
29 scale_factor = r_mean/r_mean_old
30 cube.radius *= scale_factor
31 cube.L *= scale_factor
32 cube.x *= scale_factor
33 cube.x[:,0] += numpy.abs(numpy.min(cube.x[:,0] - cube.radius))
34 cube.x[:,1] += numpy.abs(numpy.min(cube.x[:,1] - cube.radius))
35 cube.x[:,2] += numpy.abs(numpy.min(cube.x[:,2] - cube.radius))
36 cube.L[0] = numpy.max(numpy.max(cube.x[:,0] + cube.radius))
37 cube.L[1] = numpy.max(numpy.max(cube.x[:,1] + cube.radius))
38 cube.L[2] = numpy.max(numpy.max(cube.x[:,2] + cube.radius))
39 cube.writeVTK()
40
41 # Fill out grid with cubic packages
42 grid = numpy.array((
43 [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
44 [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
45 [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
46 [ 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
47 [ 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
48 [ 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
49 [ 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
50 [ 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
51 [ 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
52 [ 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
53 [ 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
54 [ 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
55 [ 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
56 [ 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
57 [ 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
58 [ 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
59 [ 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
60 [ 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
61 [ 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
62 [ 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
63 [ 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]))
64
65 # World dimensions and cube grid
66 nx = grid.shape[1] # horizontal cubes
67 ny = 1
68 nz = grid.shape[0] # vertical cubes
69 dx = cube.L[0]
70 dy = cube.L[1]
71 dz = cube.L[2]
72 Lx = dx*nx
73 Ly = dy*ny
74 Lz = dz*nz
75
76 sim = sphere.sim('cap2-cohesion=' + str(cohesion), nw=0)
77 sim.num = numpy.array([Lx, Ly, Lz])
78
79 for z in range(nz):
80 for y in range(ny):
81 for x in range(nx):
82
83 if (grid[z,x] == 0):
84 continue # skip to next iteration
85
86 for i in range(cube.np):
87 # x=x, y=y, z=z
88 pos = [ cube.x[i,0] + x*dx,
89 cube.x[i,1] + y*dy,
90 cube.x[i,2] + z*dz ]
91 sim.addParticle(pos, radius=cube.radius[i], color=grid[z,x])
92
93 cellsize_min = 2.1 * numpy.amax(sim.radius)
94 sim.defineWorldBoundaries([Lx, Ly, Lz], dx = cellsize_min)
95 sim.zeroKinematics()
96 sim.checkerboardColors()
97 sim.defaultParams(capillaryCohesion=cohesion)
98 sim.k_n[0] = 1.0e6
99 sim.k_t[0] = 1.0e6
100 sim.g[2] = -10.0
101 sim.initTemporal(2.0, epsilon=0.07)
102 sim.writeVTK()
103 sim.run(device=device)
104
105 sim.writeVTKall()