tcontactmodel_wall.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_wall.py (3777B)
---
1 #!/usr/bin/env python
2 '''
3 Validate the implemented contact models by observing the behavior of one or two
4 particles.
5 '''
6
7 import sphere
8 import numpy
9 import pytestutils
10
11 ### Wall-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=1, nw=0, sid='contactmodeltest')
17 sphere.cleanup(orig)
18 orig.radius[:] = 1.0
19 orig.x[0,:] = [5.0, 5.0, 1.05]
20 orig.vel[0,2] = -0.1
21 orig.defineWorldBoundaries(L=[10,10,10])
22 orig.gamma_wn[0] = 0.0 # Disable wall viscosity
23 orig.gamma_wt[0] = 0.0 # Disable wall viscosity
24 orig.initTemporal(total = 1.0, file_dt = 0.01)
25 #orig.time_dt = orig.time_dt*0.1
26 moment_before = orig.totalKineticEnergy()
27 orig.run(verbose=False)
28 #orig.writeVTKall()
29 orig.readlast(verbose=False)
30 pytestutils.compareFloats(orig.vel[0,2], 0.1,\
31 "Elastic normal wall collision (1/2):")
32 moment_after = orig.totalKineticEnergy()
33 #print(moment_before)
34 #print(moment_after)
35 #print("time step: " + str(orig.time_dt[0]))
36 #print(str((moment_after[0]-moment_before[0])/moment_before[0]*100.0) + " %")
37 pytestutils.compareFloats(moment_before, moment_after,\
38 "Elastic normal wall collision (2/2):")
39
40 # Oblique impact: Check for conservation of momentum (sum(v_i*m_i))
41 orig = sphere.sim(np=1, sid='contactmodeltest')
42 orig.radius[:] = 1.0
43 orig.x[0,:] = [5.0, 5.0, 1.05]
44 orig.vel[0,2] = -0.1
45 orig.vel[0,0] = 0.1
46 orig.defineWorldBoundaries(L=[10,10,10])
47 orig.gamma_wn[0] = 0.0 # Disable wall viscosity
48 orig.gamma_wt[0] = 0.0 # Disable wall viscosity
49 orig.initTemporal(total = 0.3, file_dt = 0.01)
50 moment_before = orig.totalKineticEnergy()
51 orig.run(verbose=False)
52 #orig.writeVTKall()
53 orig.readlast(verbose=False)
54 moment_after = orig.totalKineticEnergy()
55 pytestutils.compareFloats(moment_before, moment_after,\
56 " 45 deg. wall collision:\t")
57
58 ## Visco-elastic collisions
59
60 # Normal impact with normal viscous damping. Test that the lost kinetic energy
61 # is saved as dissipated viscous energy
62 orig = sphere.sim(np=1, sid='contactmodeltest')
63 orig.radius[:] = 1.0
64 orig.x[0,:] = [5.0, 5.0, 1.05]
65 orig.vel[0,2] = -0.1
66 orig.defineWorldBoundaries(L=[10,10,10])
67 orig.gamma_wn[0] = 1.0e6
68 orig.gamma_wt[0] = 0.0
69 orig.initTemporal(total = 1.0, file_dt = 0.01)
70 Ekin_before = orig.energy('kin')
71 orig.run(verbose=False)
72 #orig.writeVTKall()
73 orig.readlast(verbose=False)
74 Ekin_after = orig.energy('kin')
75 Ev_after = orig.energy('visc_n')
76 #print("Ekin_before = " + str(Ekin_before) + " J")
77 #print("Ekin_after = " + str(Ekin_after) + " J")
78 pytestutils.test(Ekin_before > Ekin_after,
79 "Viscoelastic normal wall collision (1/2):")
80 pytestutils.compareFloats(Ekin_before, Ekin_after+Ev_after,\
81 "Viscoelastic normal wall collision (2/2):", tolerance=0.05)
82
83 # Oblique impact: Check for conservation of momentum (sum(v_i*m_i))
84 orig = sphere.sim(np=1, sid='contactmodeltest')
85 orig.radius[:] = 1.0
86 orig.x[0,:] = [5.0, 5.0, 1.05]
87 orig.vel[0,2] = -0.1
88 orig.vel[0,0] = 0.1
89 orig.defineWorldBoundaries(L=[10,10,10])
90 orig.gamma_wn[0] = 1.0e6
91 orig.gamma_wt[0] = 1.0e6
92 orig.initTemporal(total = 1.0, file_dt = 0.01)
93 E_kin_before = orig.energy('kin')
94 orig.run(verbose=False)
95 #orig.writeVTKall()
96 orig.readlast(verbose=False)
97 #Ekin_after = orig.energy('kin')
98 #Erot_after = orig.energy('rot')
99 #Es_after = orig.energy('shear')
100 #pytestutils.compareFloats(Ekin_before,\
101 #Ekin_after+Erot_after+Es_after,\
102 #" 45 deg. wall collision:", tolerance=0.03)
103 pytestutils.test(Ekin_before > Ekin_after,
104 " 45 deg. wall collision (1/2):")
105 pytestutils.test((orig.angvel[0,0] == 0.0 and orig.angvel[0,1] > 0.0 \
106 and orig.angvel[0,2] == 0.0),
107 " 45 deg. wall collision (2/2):")
108
109
110
111 sphere.cleanup(orig)