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