tchannel-wet.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
---
tchannel-wet.py (6958B)
---
1 #!/usr/bin/env python
2 import sphere
3 import numpy
4
5 relaxation = False
6 consolidation = True
7 #water = False
8 water = True
9
10 id_prefix = 'chan'
11
12 #N = 2.5e3
13 #N = 5e3
14 #N = 7.5e3
15 N = 10e3
16 #N = 15e3
17 #N = 20e3
18 #N = 25e3
19 #N = 30e3
20 #N = 40e3
21
22 #dpdx = 10 # fluid-pressure gradient in Pa/m along x
23 #dpdx = 100 # fluid-pressure gradient in Pa/m along x
24 #dpdx = 200 # fluid-pressure gradient in Pa/m along x
25 dpdx = 1000
26 #dpdx = 5e3
27 #dpdx = 10e3
28 #dpdx = 20e3
29 #dpdx = 40e3
30
31 sim = sphere.sim(id_prefix + '-relax', nw=0)
32
33 if relaxation:
34 cube = sphere.sim('cube-init')
35 cube.readlast()
36 cube.adjustUpperWall(z_adjust=1.0)
37
38 # Fill out grid with cubic packages
39 grid = numpy.array((
40 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
41 [1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
42 [1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
43 [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
44 [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
45 ))
46
47 # World dimensions and cube grid
48 nx = grid.shape[1] # horizontal cubes
49 ny = 1 # horizontal (thickness) cubes
50 nz = grid.shape[0] # vertical cubes
51 dx = cube.L[0]
52 dy = cube.L[1]
53 dz = cube.L[2]
54 Lx = dx*nx
55 Ly = dy*ny
56 Lz = dz*nz
57
58 # insert particles into each cube
59 for z in range(nz):
60 for y in range(ny):
61 for x in range(nx):
62
63 if (grid[z, x] == 0):
64 continue # skip to next iteration
65
66 for i in range(cube.np):
67 pos = [cube.x[i, 0] + x*dx,
68 cube.x[i, 1] + y*dy,
69 cube.x[i, 2] + (nz-z)*dz]
70
71 sim.addParticle(pos, radius=cube.radius[i], color=grid[z, x])
72
73 # move to x=0
74 min_x = numpy.min(sim.x[:, 0] - sim.radius[:])
75 sim.x[:, 0] = sim.x[:, 0] - min_x
76
77 # move to y=0
78 min_y = numpy.min(sim.x[:, 1] - sim.radius[:])
79 sim.x[:, 1] = sim.x[:, 1] - min_y
80
81 # move to z=0
82 min_z = numpy.min(sim.x[:, 2] - sim.radius[:])
83 sim.x[:, 2] = sim.x[:, 2] - min_z
84
85 sim.defineWorldBoundaries(L=[numpy.max(sim.x[:, 0] + sim.radius[:]),
86 numpy.max(sim.x[:, 1] + sim.radius[:]),
87 numpy.max(sim.x[:, 2] + sim.radius[:])*1.2])
88 #numpy.max(sim.x[:, 2] + sim.radius[:])*10.2])
89 sim.k_t[0] = 2.0/3.0*sim.k_n[0]
90
91 # sim.cleanup()
92 sim.writeVTK()
93 print("Number of particles: " + str(sim.np))
94
95 # Set grain contact properties
96 #sim.setStiffnessNormal(1.16e7)
97 #sim.setStiffnessTangential(1.16e7)
98 sim.setYoungsModulus(70e7)
99 sim.setStaticFriction(0.5)
100 sim.setDynamicFriction(0.5)
101 sim.setDampingNormal(0.0)
102 sim.setDampingTangential(0.0)
103
104 # Set wall properties
105 sim.gamma_wn[0] = 0.0
106 sim.gamma_wt[0] = 0.0
107
108
109 # Relaxation
110
111 # Add gravitational acceleration
112 sim.g[0] = 0.0
113 sim.g[1] = 0.0
114 sim.g[2] = -9.81
115
116 sim.normalBoundariesXY()
117 # sim.consolidate(normal_stress=0.0)
118
119 # assign automatic colors, overwriting values from grid array
120 sim.checkerboardColors(nx=grid.shape[1], ny=2, nz=grid.shape[0]/4)
121
122 sim.contactmodel[0] = 2
123 sim.mu_s[0] = 0.5
124 sim.mu_d[0] = 0.5
125
126 # Set duration of simulation, automatically determine timestep, etc.
127 sim.initTemporal(total=8.0, file_dt=0.01, epsilon=0.07)
128 #sim.time_dt[0] = 1.0e-20
129 #sim.time_file_dt = sim.time_dt
130 #sim.time_total = sim.time_file_dt*5.
131 sim.zeroKinematics()
132
133 sim.run(dry=True)
134 sim.run()
135 sim.writeVTKall()
136
137
138 # Consolidation under constant normal stress
139 if consolidation:
140 sim.readlast()
141 sim.id(id_prefix + '-' + str(int(N/1000.)) + 'kPa')
142 #sim.cleanup()
143 sim.initTemporal(current=0.0, total=10.0, file_dt=0.01, epsilon=0.07)
144
145 # fix horizontal movement of lowest plane of particles
146 I = numpy.nonzero(sim.x[:, 2] < 1.5*numpy.mean(sim.radius))
147 sim.fixvel[I] = 1
148 sim.color[I] = 0
149
150 # fix horizontal movement of uppermost plane of particles
151 z_min = numpy.min(sim.x[:,2] - sim.radius)
152 z_max = numpy.max(sim.x[:,2] + sim.radius)
153 d_max_top = numpy.max(sim.radius[numpy.nonzero(sim.x[:,2] >
154 (z_max-z_min)*0.7)])*2.0
155 I = numpy.nonzero(sim.x[:,2] > (z_max - 4.0*d_max_top))
156 sim.fixvel[I] = 1
157 sim.color[I] = 0
158
159 sim.zeroKinematics()
160
161 # Wall parameters
162 sim.mu_ws[0] = 0.5
163 sim.mu_wd[0] = 0.5
164 #sim.gamma_wn[0] = 1.0e2
165 #sim.gamma_wt[0] = 1.0e2
166 sim.gamma_wn[0] = 0.0
167 sim.gamma_wt[0] = 0.0
168
169 # Particle parameters
170 sim.setYoungsModulus(70e7)
171 sim.mu_s[0] = 0.5
172 sim.mu_d[0] = 0.5
173 sim.gamma_n[0] = 0.0
174 sim.gamma_t[0] = 0.0
175
176 # apply effective normal stress from upper wall
177 sim.consolidate(normal_stress=N)
178
179
180 if water:
181
182 # read last output from previous dry experiment
183 sim.readlast()
184 sim.zeroKinematics()
185 sim.initTemporal(current=0.0, total=10.0, file_dt=0.01, epsilon=0.07)
186
187 # initialize fluid
188 sim.num = sim.num/2
189 sim.initFluid(mu=1.797e-6, p=0.0, rho=1000.0, cfd_solver=1) # water at 0 C / 1000
190 sim.setFluidCompressibility(1.426e-8) # water at 0 C
191 sim.setMaxIterations(2e5)
192 #sim.setPermeabilityGrainSize()
193 sim.setPermeabilityPrefactor(3.5e-13)
194 #sim.setPermeabilityPrefactor(3.5e-11)
195 #sim.setPermeabilityPrefactor(3.5e-15)
196
197 # initialize linear fluid pressure gradient along x
198 dx = sim.L[0]/sim.num[0]
199 for ix in numpy.arange(sim.num[0]):
200 x = dx*ix + 0.5*dx
201 sim.p_f[ix,:,:] = (dpdx*sim.L[0]) - x*dpdx
202
203 ## Fluid phase boundary conditions
204
205 # set initial pressure to zero (hydrostatic pres. distr.) everywhere
206 #sim.p_f[:,:,:] = 0.
207
208 # x
209 sim.bc_xn[0] = 0 # -x boundary: fixed pressure
210
211 # set higher fluid pressure at x-boundary away from the channel
212 #sim.p_f[0,:,:] = dpdx/sim.L[0]
213
214 sim.id(sim.id() + '-dpdx=' + str(dpdx))
215
216 sim.bc_xp[0] = 1 # +x boundary: no flow
217
218 # pressure held constant in cells at (x=nx-1, z=nz-1)
219 #sim.p_f_constant[30:,:,-2:] = 1
220 #sim.p_f_constant[40:,:,-3] = 1
221
222 # y: Can't prescribe Y pressures without affecting pressure field from
223 # -x to +x
224 #sim.setFluidYFixedPressure()
225 sim.setFluidYNoFlow()
226
227 # z
228 #sim.setFluidTopNoFlow() # ignore small contribution by IBI melting
229 sim.setFluidTopFixedPressure()
230 #sim.setFluidTopFixedFlux(specific_flux=??)
231 sim.setFluidBottomNoFlow()
232
233 # Adjust fluid grid size dynamically
234 sim.adaptiveGrid()
235 #sim.id(sim.id() + '-dpdx=' + str(dpdx) + '-newBC')
236
237
238 #sim.time_file_dt = sim.time_dt
239 #sim.time_total = sim.time_file_dt * 10
240
241 sim.run(dry=True)
242 sim.run()
243 sim.writeVTKall()