tcontactmodel.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
---
tcontactmodel.py (7660B)
---
1 #!/usr/bin/env python
2 '''
3 Validate the implemented contact models by observing the behavior of two
4 particles.
5 '''
6
7 import sphere
8 import numpy
9 import pytestutils
10
11 ### Particle-particle interaction ##############################################
12
13 ## Linear elastic collisions
14
15 # Normal impact: Check for conservation of momentum (sum(v_i*m_i))
16 orig = sphere.sim(np=2, sid='contactmodeltest')
17 after = sphere.sim(np=2, sid='contactmodeltest')
18 sphere.cleanup(orig)
19 #orig.radius[:] = [1.0, 2.0]
20 orig.radius[:] = [1.0, 1.0]
21 orig.x[0,:] = [5.0, 5.0, 2.0]
22 orig.x[1,:] = [5.0, 5.0, 4.05]
23 v_orig = 1
24 orig.vel[0,2] = v_orig
25 orig.defineWorldBoundaries(L=[10,10,10])
26 orig.initTemporal(total = 0.1, file_dt = 0.01)
27
28 orig.run(verbose=False)
29 after.readlast(verbose=False)
30 pytestutils.compareFloats(orig.vel[0,2], after.vel[1,2],\
31 "Elastic normal collision (1/4):")
32 #print(orig.totalKineticEnergy())
33 #print(after.totalKineticEnergy())
34 pytestutils.compareFloats(orig.totalKineticEnergy(), after.totalKineticEnergy(),\
35 "Elastic normal collision (2/4):")
36
37 # Normal impact with different sizes: Check for conservation of momentum
38 orig = sphere.sim(np=2, sid='contactmodeltest')
39 after = sphere.sim(np=2, sid='contactmodeltest')
40 sphere.cleanup(orig)
41 orig.radius[:] = [2.0, 1.0]
42 orig.x[0,:] = [5.0, 5.0, 2.0]
43 orig.x[1,:] = [5.0, 5.0, 5.05]
44 orig.vel[0,2] = 1.0
45 orig.defineWorldBoundaries(L=[10,10,10])
46 orig.initTemporal(total = 0.1, file_dt = 0.01)
47
48 orig.run(verbose=False)
49 after.readlast(verbose=False)
50 pytestutils.compareFloats(orig.totalKineticEnergy(), after.totalKineticEnergy(),\
51 "Elastic normal collision (3/4):")
52
53 # Normal impact with different sizes: Check for conservation of momentum
54 orig = sphere.sim(np=2, sid='contactmodeltest')
55 after = sphere.sim(np=2, sid='contactmodeltest')
56 sphere.cleanup(orig)
57 orig.radius[:] = [1.0, 2.0]
58 orig.x[0,:] = [5.0, 5.0, 2.0]
59 orig.x[1,:] = [5.0, 5.0, 5.05]
60 orig.vel[0,2] = 1.0
61 orig.defineWorldBoundaries(L=[10,10,10])
62 orig.initTemporal(total = 0.1, file_dt = 0.01)
63
64 orig.run(verbose=False)
65 after.readlast(verbose=False)
66 pytestutils.compareFloats(orig.totalKineticEnergy(), after.totalKineticEnergy(),\
67 "Elastic normal collision (4/4):")
68
69
70 ## Linear viscous-elastic collisions
71
72 # Normal impact: Check for conservation of momentum (sum(v_i*m_i))
73 orig = sphere.sim(np=2, sid='contactmodeltest')
74 after = sphere.sim(np=2, sid='contactmodeltest')
75 sphere.cleanup(orig)
76 orig.radius[:] = [1.0, 1.0]
77 orig.x[0,:] = [5.0, 5.0, 2.0]
78 orig.x[1,:] = [5.0, 5.0, 4.05]
79 v_orig = 1
80 orig.vel[0,2] = v_orig
81 orig.defineWorldBoundaries(L=[10,10,10])
82 orig.initTemporal(total = 0.1, file_dt = 0.01)
83 orig.gamma_n[0] = 1.0e6
84
85 orig.run(verbose=False)
86 after.readlast(verbose=False)
87 #print(orig.totalKineticEnergy())
88 #print(after.totalKineticEnergy())
89 #print(after.totalViscousEnergy())
90 pytestutils.test(orig.vel[0,2] > after.vel[1,2],\
91 "Viscoelastic normal collision (1/4):")
92 pytestutils.compareFloats(orig.totalKineticEnergy(),
93 after.totalKineticEnergy()
94 + after.totalViscousEnergy(),
95 "Viscoelastic normal collision (2/4):", tolerance=0.05)
96
97 # Normal impact with different sizes: Check for conservation of momentum
98 orig = sphere.sim(np=2, sid='contactmodeltest')
99 after = sphere.sim(np=2, sid='contactmodeltest')
100 sphere.cleanup(orig)
101 orig.radius[:] = [2.0, 1.0]
102 orig.x[0,:] = [5.0, 5.0, 2.0]
103 orig.x[1,:] = [5.0, 5.0, 5.05]
104 orig.vel[0,2] = 1.0
105 orig.defineWorldBoundaries(L=[10,10,10])
106 orig.initTemporal(total = 0.1, file_dt = 0.01)
107 orig.gamma_n[0] = 1.0e6
108
109 orig.run(verbose=False)
110 after.readlast(verbose=False)
111 pytestutils.compareFloats(orig.totalKineticEnergy(),
112 after.totalKineticEnergy()
113 + after.totalViscousEnergy(),
114 "Viscoelastic normal collision (3/4):", tolerance=0.05)
115
116 # Normal impact with different sizes: Check for conservation of momentum
117 orig = sphere.sim(np=2, sid='contactmodeltest')
118 after = sphere.sim(np=2, sid='contactmodeltest')
119 sphere.cleanup(orig)
120 orig.radius[:] = [1.0, 2.0]
121 orig.x[0,:] = [5.0, 5.0, 2.0]
122 orig.x[1,:] = [5.0, 5.0, 5.05]
123 orig.vel[0,2] = 1.0
124 orig.defineWorldBoundaries(L=[10,10,10])
125 orig.initTemporal(total = 0.1, file_dt = 0.01)
126 orig.gamma_n[0] = 1.0e6
127
128 orig.run(verbose=False)
129 after.readlast(verbose=False)
130 pytestutils.compareFloats(orig.totalKineticEnergy(),
131 after.totalKineticEnergy()
132 + after.totalViscousEnergy(),
133 "Viscoelastic normal collision (4/4):", tolerance=0.05)
134
135
136
137 ## Oblique elastic collisions
138
139 # Normal impact, low angle, no slip
140 orig = sphere.sim(np=2, sid='contactmodeltest')
141 after = sphere.sim(np=2, sid='contactmodeltest')
142 sphere.cleanup(orig)
143 #orig.radius[:] = [1.0, 2.0]
144 orig.radius[:] = [1.0, 1.0]
145 orig.x[0,:] = [5.0, 5.0, 2.0]
146 orig.x[1,:] = [5.0, 5.0, 4.05]
147 orig.vel[0,2] = 1
148 orig.vel[0,0] = 1
149 orig.mu_s[0] = 1e9 # no slip
150 orig.mu_d[0] = 1e9 # no slip
151 orig.defineWorldBoundaries(L=[10,10,10])
152 orig.initTemporal(total = 0.1, file_dt = 0.01)
153
154 orig.run(verbose=False)
155 after.readlast(verbose=False)
156 pytestutils.test((after.angvel[:,1] < 0.0).all(),
157 "Oblique normal collision (1/8):")
158 pytestutils.compareFloats(orig.totalKineticEnergy(),
159 after.totalKineticEnergy()
160 + after.totalRotationalEnergy(),
161 "Oblique normal collision (2/8):", tolerance=0.05)
162
163 # Normal impact, low angle, slip
164 orig = sphere.sim(np=2, sid='contactmodeltest')
165 after = sphere.sim(np=2, sid='contactmodeltest')
166 sphere.cleanup(orig)
167 #orig.radius[:] = [1.0, 2.0]
168 orig.radius[:] = [1.0, 1.0]
169 orig.x[0,:] = [5.0, 5.0, 2.0]
170 orig.x[1,:] = [5.0, 5.0, 4.05]
171 orig.vel[0,2] = 1
172 orig.vel[0,0] = 1
173 orig.mu_s[0] = 0.3
174 orig.mu_d[0] = 0.3
175 orig.defineWorldBoundaries(L=[10,10,10])
176 orig.initTemporal(total = 0.1, file_dt = 0.01)
177
178 orig.run(verbose=False)
179 after.readlast(verbose=False)
180 pytestutils.compareFloats(orig.totalKineticEnergy(),
181 after.totalKineticEnergy()
182 + after.totalRotationalEnergy()
183 + after.totalFrictionalEnergy(),
184 "Oblique normal collision (3/8):", tolerance=0.05)
185 pytestutils.test((after.angvel[:,1] < 0.0).all(),
186 "Oblique normal collision (4/8):")
187 pytestutils.test(after.totalFrictionalEnergy() > 0.0,
188 "Oblique normal collision (5/8):")
189
190 # Normal impact, low angle, slip, viscous damping tangentially
191 orig = sphere.sim(np=2, sid='contactmodeltest')
192 after = sphere.sim(np=2, sid='contactmodeltest')
193 sphere.cleanup(orig)
194 #orig.radius[:] = [1.0, 2.0]
195 orig.radius[:] = [1.0, 1.0]
196 orig.x[0,:] = [5.0, 5.0, 2.0]
197 orig.x[1,:] = [5.0, 5.0, 4.05]
198 orig.vel[0,2] = 1
199 orig.vel[0,0] = 1
200 orig.mu_s[0] = 0.3
201 orig.mu_d[0] = 0.3
202 orig.gamma_t[0] = 1.0e3
203 orig.defineWorldBoundaries(L=[10,10,10])
204 orig.initTemporal(total = 0.1, file_dt = 0.01)
205
206 orig.run(verbose=False)
207 after.readlast(verbose=False)
208 print(after.totalViscousEnergy())
209 pytestutils.compareFloats(orig.totalKineticEnergy(),
210 after.totalKineticEnergy()
211 + after.totalRotationalEnergy()
212 + after.totalFrictionalEnergy()
213 + after.totalViscousEnergy(),
214 "Oblique normal collision (6/8):", tolerance=0.05)
215 pytestutils.test((after.angvel[:,1] < 0.0).all(),
216 "Oblique normal collision (7/8):")
217 pytestutils.test(after.totalFrictionalEnergy() > 0.0,
218 "Oblique normal collision (8/8):")
219 pytestutils.test(after.totalFrictionalEnergy() > 0.0,
220 "Oblique normal collision (8/8):")
221
222 orig.cleanup()