tcreep-master.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
---
tcreep-master.py (5380B)
---
1 #!/usr/bin/env python
2
3 # Import sphere functionality
4 import sphere
5 import numpy
6
7 ### EXPERIMENT SETUP ###
8 initialization = False
9 consolidation = False
10 shearing = False
11 creeping = True
12 rendering = False
13 plots = True
14
15 # Common simulation id
16 sim_id = "creep2"
17
18 # Fluid-pressure gradient [Pa/m]
19 dpdx = -100.0
20
21 # Deviatoric stress [Pa]
22 N = 100e3
23
24 # Grain density
25 rho_g = 1000.0
26
27 # Fluid density
28 rho_f = 1000.0
29
30 # Gravitational acceleration
31 g = 10.0
32
33 # Number of particles
34 np = 1e4
35
36
37 ### INITIALIZATION ###
38
39 # New class
40 #init = sphere.sim(np = np, nd = 3, nw = 0, sid = sim_id + "-init")
41 init = sphere.sim(np = np, nd = 3, nw = 0, sid = 'creep1' + "-init")
42
43 # Uniform radii from 0.8 cm to 1.2 cm
44 init.generateRadii(psd = 'uni', mean = 0.005, variance = 0.001)
45
46 # Use default params
47 init.defaultParams(gamma_n = 100.0, mu_s = 0.6, mu_d = 0.6)
48 init.setYoungsModulus(1e8)
49
50 # Add gravity
51 init.g[2] = -g
52
53 # Periodic x and y boundaries
54 init.periodicBoundariesXY()
55
56 # Initialize positions in random grid (also sets world size)
57 hcells = np**(1.0/3.0)
58 init.initRandomGridPos(gridnum = [hcells, hcells, 1e9])
59
60 # Set duration of simulation
61 init.initTemporal(total = 5.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')
72
73 init.writeVTKall()
74
75 if (rendering == True):
76 # Render images with raytracer
77 init.render(method = "angvel", max_val = 0.3, verbose = False)
78
79
80 ### CONSOLIDATION ###
81
82 # New class
83 cons = sphere.sim(np = init.np, nw = 1, sid = sim_id +
84 "-cons-N{}".format(N))
85
86 # Read last output file of initialization step
87 lastf = init.status()
88 #cons.readbin("../output/" + sim_id + "-init.output{:0=5}.bin".format(lastf), verbose=False)
89 cons.readbin("../output/" + 'creep1' + "-init.output{:0=5}.bin".format(lastf),
90 verbose=False)
91
92 # Periodic x and y boundaries
93 cons.periodicBoundariesXY()
94
95 # Setup consolidation experiment
96 cons.consolidate(normal_stress = N)
97
98 # Disable wall viscosities
99 cons.gamma_wn[0] = 0.0
100 cons.gamma_wt[0] = 0.0
101
102 cons.rho[0] = rho_g
103 cons.g[2] = -g
104
105 # Set duration of simulation
106 cons.initTemporal(total = 1.5)
107
108 if (consolidation == True):
109
110 # Run sphere
111 cons.run(dry = True) # show values, don't run
112 cons.run() # run
113
114 if (plots == True):
115 # Make a graph of energies
116 cons.visualize('energy')
117 cons.visualize('walls')
118
119 cons.writeVTKall()
120
121 if (rendering == True):
122 # Render images with raytracer
123 cons.render(method = "pres", max_val = 2.0*N, verbose = False)
124
125
126 ### SHEARING ###
127
128 # New class
129 shear = sphere.sim(np = cons.np, nw = cons.nw, sid = sim_id +
130 "-shear-N{}".format(N))
131
132 # Read last output file of initialization step
133 lastf = cons.status()
134 shear.readbin("../output/" + sim_id +
135 "-cons-N{}.output{:0=5}.bin".format(N, lastf), verbose =
136 False)
137
138 # Periodic x and y boundaries
139 shear.periodicBoundariesXY()
140
141 shear.rho[0] = rho_g
142 shear.g[2] = -g
143
144 # Disable particle viscosities
145 shear.gamma_n[0] = 0.0
146 shear.gamma_t[0] = 0.0
147
148 # Setup shear experiment
149 shear.shear(shear_strain_rate = 0.1)
150
151 # Set duration of simulation
152 shear.initTemporal(total = 10.0)
153
154 if (shearing == True):
155
156 # Run sphere
157 shear.run(dry = True)
158 shear.run()
159
160 if (plots == True):
161 # Make a graph of energies
162 shear.visualize('energy')
163 shear.visualize('shear')
164
165 shear.writeVTKall()
166
167 if (rendering == True):
168 # Render images with raytracer
169 shear.render(method = "pres", max_val = 2.0*N, verbose = False)
170
171
172 ### CREEP ###
173
174 # New class
175 creep = sphere.sim(np = cons.np, nw = cons.nw, sid = sim_id +
176 "-N{}-dpdx{}".format(N, dpdx))
177
178 # Read last output file of initialization step
179 lastf = shear.status()
180 creep.readbin("../output/" + sim_id +
181 "-shear-N{}.output{:0=5}.bin".format(N, lastf), verbose =
182 False)
183
184 # Periodic x and y boundaries
185 creep.periodicBoundariesXY()
186
187 # set fluid and solver properties
188 creep.initFluid(mu=8.9e-4, p=0.0, rho=rho_f, cfd_solver=1) # water at 25 C
189 creep.setMaxIterations(2e5)
190 creep.setPermeabilityGrainSize()
191 creep.setFluidCompressibility(4.6e-10) # water at 25 C
192
193 # set fluid BCs
194 creep.setFluidTopNoFlow()
195 creep.setFluidBottomNoFlow()
196 creep.setFluidXFixedPressure()
197 creep.adaptiveGrid()
198
199 # set fluid pressures at the boundaries and internally
200 dx = creep.L[0]/creep.num[0]
201 for ix in range(creep.num[0]):
202 x = ix + 0.5*dx
203 creep.p_f[ix,:,:] = numpy.abs(creep.L[0]*dpdx) + x*dpdx
204
205 creep.zeroKinematics()
206
207 # Remove fixvel constraint from uppermost grains
208 creep.fixvel[numpy.nonzero(creep.x[:,2] > 0.5*creep.L[2])] = 0
209
210 # Produce regular coloring pattern
211 creep.checkerboardColors(creep.num[0], creep.num[1], creep.num[2])
212 creep.color[numpy.nonzero(creep.fixvel == 1)] == -1
213
214 # Adapt grid size during progressive deformation
215 creep.adaptiveGrid()
216
217 # Set duration of simulation
218 creep.initTemporal(total = 20.0)
219
220 if (creeping == True):
221
222 # Run sphere
223 creep.run(dry = True)
224 creep.run()
225
226 if (plots == True):
227 # Make a graph of energies
228 creep.visualize('energy')
229
230 creep.writeVTKall()
231
232 if (rendering == True):
233 # Render images with raytracer
234 creep.render(method = "pres", max_val = 2.0*N, verbose = False)