 |
tchannel-wet.py - sphere - GPU-based 3D discrete element method algorithm with … |
 |
git clone git://src.adamsgaard.dk/sphere |
 |
Log |
|
Files |
|
Refs |
|
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=gr… |
|
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… |
|
88 #numpy.max(sim.x[:, 2] + sim.radius[:])*… |
|
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)])… |
|
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=… |
|
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) # wa… |
|
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.) everyw… |
|
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 fiel… |
|
223 # -x to +x |
|
224 #sim.setFluidYFixedPressure() |
|
225 sim.setFluidYNoFlow() |
|
226 |
|
227 # z |
|
228 #sim.setFluidTopNoFlow() # ignore small contribution by IBI mel… |
|
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() |
