tcohesion.jl - Granular.jl - Julia package for granular dynamics simulation
(HTM) git clone git://src.adamsgaard.dk/Granular.jl
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) LICENSE
---
tcohesion.jl (8100B)
---
1 #!/usr/bin/env julia
2 using Test
3 import Granular
4
5 # Check for conservation of kinetic energy (=momentum) during a normal collision
6 # between two ice cylindrical grains
7
8 verbose=false
9
10 sim_init = Granular.createSimulation()
11 Granular.addGrainCylindrical!(sim_init, [0., 0.], 10., 1.)
12 Granular.addGrainCylindrical!(sim_init, [18., 0.], 10., 1.)
13 sim_init.grains[1].youngs_modulus = 1e-5 # repulsion is negligible
14 sim_init.grains[2].youngs_modulus = 1e-5 # repulsion is negligible
15 Granular.setTimeStep!(sim_init, verbose=verbose)
16
17 @info "# Check contact age scheme"
18 sim = deepcopy(sim_init)
19 Granular.setTotalTime!(sim, 10.)
20 sim.time_step = 1.
21 Granular.run!(sim, verbose=verbose)
22 Granular.removeSimulationFiles(sim)
23 @test sim.grains[1].contact_age[1] ≈ sim.time
24
25 @info "# Check if bonds add tensile strength"
26 sim = Granular.createSimulation(id="cohesion")
27 Granular.addGrainCylindrical!(sim, [0., 0.], 10., 1., tensile_strength=500e3)
28 Granular.addGrainCylindrical!(sim, [20.1, 0.], 10., 1., tensile_strength=500e3)
29 sim.grains[1].lin_vel[1] = 0.1
30 Granular.setTimeStep!(sim)
31 Granular.setTotalTime!(sim, 10.)
32 Granular.run!(sim, verbose=verbose)
33 Granular.removeSimulationFiles(sim)
34 @test sim.grains[1].lin_vel[1] > 0.
35 @test sim.grains[1].lin_vel[2] ≈ 0.
36 @test sim.grains[2].lin_vel[1] > 0.
37 @test sim.grains[2].lin_vel[2] ≈ 0.
38 @test sim.grains[1].ang_vel ≈ zeros(3)
39 @test sim.grains[2].ang_vel ≈ zeros(3)
40
41 @info "# Add shear strength and test bending resistance (one grain rotating)"
42 sim = Granular.createSimulation(id="cohesion")
43 Granular.addGrainCylindrical!(sim, [0., 0.], 10.1, 1., tensile_strength=500e3,
44 shear_strength=500e3)
45 Granular.addGrainCylindrical!(sim, [20., 0.], 10., 1., tensile_strength=500e3,
46 shear_strength=500e3)
47 sim.grains[1].ang_vel[3] = 0.1
48 Granular.findContacts!(sim) # make sure contact is registered
49 sim.grains[1].contact_radius=10.0 # decrease radius so there isn't compression
50 E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim)
51 E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim)
52 Granular.setTimeStep!(sim)
53 Granular.setTotalTime!(sim, 5.)
54 #Granular.setOutputFileInterval!(sim, 0.1)
55 Granular.run!(sim, verbose=verbose)
56 Granular.removeSimulationFiles(sim)
57 @test sim.grains[1].lin_vel[1] ≈ 0.
58 @test sim.grains[1].lin_vel[2] ≈ 0.
59 @test sim.grains[2].lin_vel[1] ≈ 0.
60 @test sim.grains[2].lin_vel[2] ≈ 0.
61 @test sim.grains[1].ang_vel[3] != 0.
62 @test sim.grains[2].ang_vel[3] != 0.
63 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
64 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
65 E_therm_final = Granular.totalGrainThermalEnergy(sim)
66 @test E_kin_lin_init ≈ E_kin_lin_final
67 @test E_kin_rot_init > E_kin_rot_final + E_therm_final
68
69 @info "# Add shear strength and test bending resistance (one grain rotating)"
70 sim = Granular.createSimulation(id="cohesion")
71 Granular.addGrainCylindrical!(sim, [0., 0.], 10.1, 1., tensile_strength=500e3,
72 shear_strength=500e3)
73 Granular.addGrainCylindrical!(sim, [20., 0.], 10., 1., tensile_strength=500e3,
74 shear_strength=500e3)
75 sim.grains[2].ang_vel[3] = 0.1
76 Granular.findContacts!(sim) # make sure contact is registered
77 sim.grains[1].contact_radius=10.0 # decrease radius so there isn't compression
78 E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim)
79 E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim)
80 Granular.setTimeStep!(sim)
81 Granular.setTotalTime!(sim, 5.)
82 #Granular.setOutputFileInterval!(sim, 0.1)
83 Granular.run!(sim, verbose=verbose)
84 Granular.removeSimulationFiles(sim)
85 @test sim.grains[1].lin_vel[1] ≈ 0.
86 @test sim.grains[1].lin_vel[2] ≈ 0.
87 @test sim.grains[2].lin_vel[1] ≈ 0.
88 @test sim.grains[2].lin_vel[2] ≈ 0.
89 @test sim.grains[1].ang_vel[3] != 0.
90 @test sim.grains[2].ang_vel[3] != 0.
91 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
92 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
93 E_therm_final = Granular.totalGrainThermalEnergy(sim)
94 @test E_kin_lin_init ≈ E_kin_lin_final
95 @test E_kin_rot_init > E_kin_rot_final + E_therm_final
96
97 @info "# Add shear strength and test bending resistance (both grains rotating)"
98 sim = Granular.createSimulation(id="cohesion")
99 Granular.addGrainCylindrical!(sim, [0., 0.], 10.0000001, 1., tensile_strength=500e3,
100 shear_strength=500e3)
101 Granular.addGrainCylindrical!(sim, [20., 0.], 10., 1., tensile_strength=500e3,
102 shear_strength=500e3)
103 sim.grains[1].ang_vel[3] = 0.1
104 sim.grains[2].ang_vel[3] = -0.1
105 Granular.findContacts!(sim) # make sure contact is registered
106 sim.grains[1].contact_radius=10.0 # decrease radius so there isn't compression
107 E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim)
108 E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim)
109 Granular.setTimeStep!(sim)
110 Granular.setTotalTime!(sim, 5.)
111 #Granular.setOutputFileInterval!(sim, 0.1)
112 Granular.run!(sim, verbose=verbose)
113 Granular.removeSimulationFiles(sim)
114 @test sim.grains[1].lin_vel[1] ≈ 0.
115 @test sim.grains[1].lin_vel[2] ≈ 0.
116 @test sim.grains[2].lin_vel[1] ≈ 0.
117 @test sim.grains[2].lin_vel[2] ≈ 0.
118 @test sim.grains[1].ang_vel[3] != 0.
119 @test sim.grains[2].ang_vel[3] != 0.
120 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
121 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
122 E_therm_final = Granular.totalGrainThermalEnergy(sim)
123 @test E_kin_lin_init ≈ E_kin_lin_final
124 @test E_kin_rot_init > E_kin_rot_final + E_therm_final
125
126 @info "# Break bond through bending I"
127 sim = Granular.createSimulation(id="cohesion")
128 Granular.addGrainCylindrical!(sim, [0., 0.], 10.0000001, 1., tensile_strength=500e3,
129 shear_strength=500e3)
130 Granular.addGrainCylindrical!(sim, [20., 0.], 10., 1., tensile_strength=500e3,
131 shear_strength=500e3)
132 sim.grains[1].ang_vel[3] = 100
133 sim.grains[2].ang_vel[3] = -100
134 Granular.findContacts!(sim) # make sure contact is registered
135 sim.grains[1].contact_radius=10.0 # decrease radius so there isn't compression
136 E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim)
137 E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim)
138 Granular.setTimeStep!(sim)
139 Granular.setTotalTime!(sim, 5.)
140 #Granular.setOutputFileInterval!(sim, 0.1)
141 Granular.run!(sim, verbose=verbose)
142 Granular.removeSimulationFiles(sim)
143 @test sim.grains[1].lin_vel[1] ≈ 0.
144 @test sim.grains[1].lin_vel[2] ≈ 0.
145 @test sim.grains[2].lin_vel[1] ≈ 0.
146 @test sim.grains[2].lin_vel[2] ≈ 0.
147 @test sim.grains[1].ang_vel[3] != 0.
148 @test sim.grains[2].ang_vel[3] != 0.
149 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
150 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
151 E_therm_final = Granular.totalGrainThermalEnergy(sim)
152 @test E_kin_lin_init ≈ E_kin_lin_final
153 @test sim.grains[1].n_contacts == 0
154 @test sim.grains[2].n_contacts == 0
155
156 @info "# Break bond through bending II"
157 sim = Granular.createSimulation(id="cohesion")
158 Granular.addGrainCylindrical!(sim, [0., 0.], 10.1, 1., tensile_strength=500e3,
159 shear_strength=50e3)
160 Granular.addGrainCylindrical!(sim, [20., 0.], 10., 1., tensile_strength=500e3,
161 shear_strength=50e3)
162 sim.grains[1].ang_vel[3] = 100
163 sim.grains[2].ang_vel[3] = 0.0
164 Granular.findContacts!(sim) # make sure contact is registered
165 sim.grains[1].contact_radius=10.0 # decrease radius so there isn't compression
166 E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim)
167 E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim)
168 Granular.setTimeStep!(sim)
169 Granular.setTotalTime!(sim, 5.)
170 #Granular.setOutputFileInterval!(sim, 0.1)
171 Granular.run!(sim, verbose=verbose)
172 Granular.removeSimulationFiles(sim)
173 @test sim.grains[1].lin_vel[1] ≈ 0.
174 @test sim.grains[1].lin_vel[2] ≈ 0.
175 @test sim.grains[2].lin_vel[1] ≈ 0.
176 @test sim.grains[2].lin_vel[2] ≈ 0.
177 @test sim.grains[1].ang_vel[3] != 0.
178 @test sim.grains[2].ang_vel[3] != 0.
179 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
180 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
181 E_therm_final = Granular.totalGrainThermalEnergy(sim)
182 @test E_kin_lin_init ≈ E_kin_lin_final
183 @test sim.grains[1].n_contacts == 0
184 @test sim.grains[2].n_contacts == 0