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