tcfd_tests_darcy_particles.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
---
tcfd_tests_darcy_particles.py (6788B)
---
1 #!/usr/bin/env python
2 from pytestutils import *
3 import sphere
4 import numpy
5
6 #'''
7 print("### Steady state, no gravity, no forcing, Dirichlet+Dirichlet BCs")
8 orig = sphere.sim('darcy_particles', np = 1000)
9 orig.cleanup()
10 #orig.generateRadii(histogram=False, psd='uni', radius_mean=5.0e-4, radius_variance=5.0e-5)
11 orig.defaultParams()
12 orig.generateRadii(psd='uni', mean=5.0e-2, variance=5.0e-5)
13 orig.initRandomGridPos([20, 20, 200])
14 orig.initTemporal(total=0.005, file_dt=0.001)
15 orig.initFluid(cfd_solver=1)
16 #orig.p_f[5,3,2] *= 1.5
17 #orig.k_c[0] = 4.6e-15
18 orig.k_c[0] = 4.6e-10
19 #orig.g[2] = -10.0
20 orig.setStiffnessNormal(36.4e9)
21 orig.setStiffnessTangential(36.4e9/3.0)
22 orig.run(verbose=False)
23 #orig.writeVTKall()
24 py = sphere.sim(sid = orig.sid, fluid = True)
25 py.readlast(verbose=False)
26
27 zeros = numpy.zeros((orig.num))
28 py.readlast(verbose = False)
29 compareNumpyArrays(zeros, py.p_f, "Conservation of pressure:")
30
31 # Fluid flow should be very small
32 if ((numpy.abs(py.v_f[:,:,:,:]) < 1.0e-6).all()):
33 print("Flow field:\t\t" + passed())
34 else:
35 print("Flow field:\t\t" + failed())
36 print(numpy.min(py.v_f))
37 print(numpy.mean(py.v_f))
38 print(numpy.max(py.v_f))
39 raise Exception("Failed")
40
41
42
43 print("### Steady state, no gravity, no forcing, Neumann+Dirichlet BCs")
44 orig = sphere.sim('darcy_particles', np = 1000)
45 orig.cleanup()
46 #orig.generateRadii(histogram=False, psd='uni', radius_mean=5.0e-4, radius_variance=5.0e-5)
47 orig.defaultParams()
48 orig.generateRadii(psd='uni', mean=5.0e-2, variance=5.0e-5)
49 orig.initRandomGridPos([20, 20, 200])
50 orig.initTemporal(total=0.005, file_dt=0.001)
51 orig.initFluid(cfd_solver=1)
52 #orig.p_f[5,3,2] *= 1.5
53 #orig.k_c[0] = 4.6e-15
54 orig.k_c[0] = 4.6e-10
55 orig.setFluidBottomNoFlow()
56 #orig.g[2] = -10.0
57 orig.setStiffnessNormal(36.4e9)
58 orig.setStiffnessTangential(36.4e9/3.0)
59 orig.run(verbose=False)
60 #orig.writeVTKall()
61 py = sphere.sim(sid = orig.sid, fluid = True)
62 py.readlast(verbose=False)
63
64 zeros = numpy.zeros((orig.num))
65 py.readlast(verbose = False)
66 compareNumpyArrays(zeros, py.p_f, "Conservation of pressure:")
67
68 # Fluid flow should be very small
69 if ((numpy.abs(py.v_f[:,:,:,:]) < 1.0e-6).all()):
70 print("Flow field:\t\t" + passed())
71 else:
72 print("Flow field:\t\t" + failed())
73 print(numpy.min(py.v_f))
74 print(numpy.mean(py.v_f))
75 print(numpy.max(py.v_f))
76 raise Exception("Failed")
77
78
79
80 print("### Steady state, no gravity, no forcing, Neumann+Neumann BCs")
81 orig = sphere.sim('darcy_particles', np = 1000)
82 orig.cleanup()
83 #orig.generateRadii(histogram=False, psd='uni', radius_mean=5.0e-4, radius_variance=5.0e-5)
84 orig.defaultParams()
85 orig.generateRadii(psd='uni', mean=5.0e-2, variance=5.0e-5)
86 orig.initRandomGridPos([20, 20, 200])
87 orig.initTemporal(total=0.005, file_dt=0.001)
88 orig.initFluid(cfd_solver=1)
89 #orig.p_f[5,3,2] *= 1.5
90 #orig.k_c[0] = 4.6e-15
91 orig.k_c[0] = 4.6e-10
92 orig.setFluidBottomNoFlow()
93 orig.setFluidTopNoFlow()
94 #orig.g[2] = -10.0
95 orig.setStiffnessNormal(36.4e9)
96 orig.setStiffnessTangential(36.4e9/3.0)
97 orig.run(verbose=False)
98 #orig.writeVTKall()
99 py = sphere.sim(sid = orig.sid, fluid = True)
100 py.readlast(verbose=False)
101
102 zeros = numpy.zeros((orig.num))
103 py.readlast(verbose = False)
104 compareNumpyArrays(zeros, py.p_f, "Conservation of pressure:")
105
106 # Fluid flow should be very small
107 if ((numpy.abs(py.v_f[:,:,:,:]) < 1.0e-6).all()):
108 print("Flow field:\t\t" + passed())
109 else:
110 print("Flow field:\t\t" + failed())
111 print(numpy.min(py.v_f))
112 print(numpy.mean(py.v_f))
113 print(numpy.max(py.v_f))
114 raise Exception("Failed")
115 #'''
116
117
118 print("### Dynamic wall: Transient, gravity, Dirichlet+Neumann BCs")
119 #orig = sphere.sim('diffusivity-relax', fluid=False)
120 orig = sphere.sim('cube-init', fluid=False)
121 orig.readlast(verbose=False)
122 orig.num[2] /= 2
123 orig.L[2] /= 2.0
124 orig.id('darcy_fluidization')
125 orig.cleanup()
126 orig.setStiffnessNormal(36.4e9)
127 orig.setStiffnessTangential(36.4e9/3.0)
128 orig.initTemporal(total=0.0005, file_dt=0.0001)
129 orig.initFluid(cfd_solver=1)
130 orig.setFluidBottomNoFlow()
131 orig.g[2] = -10.0
132 #orig.k_c[0] = numpy.mean(orig.radius)**2/540.0
133 orig.consolidate()
134
135 orig.run(verbose=False)
136 #orig.writeVTKall()
137 py = sphere.sim(sid = orig.sid, fluid = True)
138 py.readlast(verbose=False)
139 test(orig.w_x[0] > py.w_x[0], 'Wall movement:\t\t')
140
141 print("### Dynamic wall: Transient, gravity, Dirichlet+Neumann BCs, ndem=10")
142 #orig = sphere.sim('diffusivity-relax', fluid=False)
143 orig = sphere.sim('cube-init', fluid=False)
144 orig.readlast(verbose=False)
145 orig.num[2] /= 2
146 orig.L[2] /= 2.0
147 orig.id('darcy_fluidization')
148 orig.cleanup()
149 orig.setStiffnessNormal(36.4e9)
150 orig.setStiffnessTangential(36.4e9/3.0)
151 orig.initTemporal(total=0.0005, file_dt=0.0001)
152 orig.initFluid(cfd_solver=1)
153 orig.setDEMstepsPerCFDstep(10)
154 orig.setFluidBottomNoFlow()
155 orig.g[2] = -10.0
156 #orig.k_c[0] = numpy.mean(orig.radius)**2/540.0
157 orig.consolidate()
158
159 orig.run(verbose=False)
160 #orig.writeVTKall()
161 py = sphere.sim(sid = orig.sid, fluid = True)
162 py.readlast(verbose=False)
163 test(orig.w_x[0] > py.w_x[0], 'Wall movement:\t\t')
164
165
166 print("### Fluidization test: Transient, gravity, Dirichlet+Dirichlet BCs")
167 #orig = sphere.sim('diffusivity-relax', fluid=False)
168 orig = sphere.sim('cube-init', fluid=False)
169 orig.readlast(verbose=False)
170 orig.num[2] /= 2
171 orig.L[2] /= 2.0
172 orig.id('darcy_fluidization')
173 orig.cleanup()
174 orig.setDEMstepsPerCFDstep(1)
175 orig.setStiffnessNormal(36.4e9)
176 orig.setStiffnessTangential(36.4e9/3.0)
177 orig.initTemporal(total=0.0005, file_dt=0.0001)
178 orig.initFluid(cfd_solver=1)
179 orig.g[2] = -10.0
180 #orig.k_c[0] = numpy.mean(orig.radius)**2/540.0
181
182 mean_porosity = orig.bulkPorosity()
183 fluidize_pressure = numpy.abs((orig.rho - orig.rho_f) \
184 *(1.0 - mean_porosity)*numpy.abs(orig.g[2]))
185
186 fluid_pressure_gradient = numpy.array([0.9, 2.0])
187
188 for i in numpy.arange(fluid_pressure_gradient.size):
189
190 orig.id('fluidization-' + str(fluid_pressure_gradient[i]))
191 # set pressure gradient
192 dpdz = fluid_pressure_gradient[i] * fluidize_pressure
193 dp = dpdz * orig.L[2]
194 base_p = 0.0
195 orig.p_f[:,:,0] = base_p + dp # high pressure at bottom
196 orig.p_f[:,:,-1] = base_p # low pressure at top
197
198 orig.run(verbose=False)
199 #orig.writeVTKall()
200 py = sphere.sim(sid = orig.sid, fluid = True)
201 py.readlast(verbose=False)
202
203 """print('Mean particle velocity: '
204 + str(numpy.mean(py.vel[:,0])) + ', '
205 + str(numpy.mean(py.vel[:,1])) + ', '
206 + str(numpy.mean(py.vel[:,2])) + ' m/s')"""
207
208 z_vel_threshold = 0.002
209 if fluid_pressure_gradient[i] < 1.0:
210 test(numpy.mean(py.vel[:,2]) < z_vel_threshold,
211 'Fluidization (' + str(fluid_pressure_gradient[i]) + '):\t')
212 elif fluid_pressure_gradient[i] > 1.0:
213 test(numpy.mean(py.vel[:,2]) > z_vel_threshold,
214 'Fluidization (' + str(fluid_pressure_gradient[i]) + '):\t')
215 orig.cleanup()