collision-2floes-normal.jl (13115B)
1 #!/usr/bin/env julia 2 3 # Check for conservation of kinetic energy (=momentum) during a normal collision 4 # between two ice cylindrical grains 5 6 verbose=false 7 8 @info "# One ice floe fixed" 9 sim = Granular.createSimulation(id="test") 10 Granular.addGrainCylindrical!(sim, [0., 0.], 10., 1., verbose=verbose) 11 Granular.addGrainCylindrical!(sim, [20.05, 0.], 10., 1., verbose=verbose) 12 sim.grains[1].lin_vel[1] = 0.1 13 sim.grains[2].fixed = true 14 15 E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim) 16 E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim) 17 18 # With decreasing timestep (epsilon towards 0), the explicit integration scheme 19 # should become more correct 20 21 Granular.setTotalTime!(sim, 10.0) 22 sim_init = deepcopy(sim) 23 24 @info "Testing kinetic energy conservation with Two-term Taylor scheme" 25 Granular.setTimeStep!(sim, epsilon=0.07) 26 tol = 0.2 27 @info "Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)" 28 Granular.run!(sim, temporal_integration_method="Two-term Taylor", verbose=verbose) 29 30 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim) 31 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim) 32 @test E_kin_lin_init ≈ E_kin_lin_final atol=E_kin_lin_init*tol 33 @test E_kin_rot_init ≈ E_kin_rot_final 34 35 36 @info "Testing kinetic energy conservation with Two-term Taylor scheme" 37 sim = deepcopy(sim_init) 38 Granular.setTimeStep!(sim, epsilon=0.007) 39 tol = 0.02 40 @info "Relative tolerance: $(tol*100.)%" 41 Granular.run!(sim, temporal_integration_method="Two-term Taylor", verbose=verbose) 42 43 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim) 44 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim) 45 @test E_kin_lin_init ≈ E_kin_lin_final atol=E_kin_lin_init*tol 46 @test E_kin_rot_init ≈ E_kin_rot_final 47 48 49 @info "Testing kinetic energy conservation with Three-term Taylor scheme" 50 sim = deepcopy(sim_init) 51 Granular.setTimeStep!(sim, epsilon=0.07) 52 tol = 0.01 53 @info "Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)" 54 Granular.run!(sim, temporal_integration_method="Three-term Taylor", 55 verbose=verbose) 56 57 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim) 58 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim) 59 @test E_kin_lin_init ≈ E_kin_lin_final atol=E_kin_lin_init*tol 60 @test E_kin_rot_init ≈ E_kin_rot_final 61 62 63 @info "# Ice floes free to move" 64 65 sim = Granular.createSimulation(id="test") 66 Granular.addGrainCylindrical!(sim, [0., 0.], 10., 1., verbose=verbose) 67 Granular.addGrainCylindrical!(sim, [20.05, 0.], 10., 1., verbose=verbose) 68 sim.grains[1].lin_vel[1] = 0.1 69 70 E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim) 71 E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim) 72 73 # With decreasing timestep (epsilon towards 0), the explicit integration scheme 74 # should become more correct 75 76 Granular.setTotalTime!(sim, 10.0) 77 sim_init = deepcopy(sim) 78 79 @info "Testing kinetic energy conservation with Two-term Taylor scheme" 80 Granular.setTimeStep!(sim, epsilon=0.07) 81 tol = 0.2 82 @info "Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)" 83 Granular.run!(sim, temporal_integration_method="Two-term Taylor", verbose=verbose) 84 85 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim) 86 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim) 87 @test E_kin_lin_init ≈ E_kin_lin_final atol=E_kin_lin_init*tol 88 @test E_kin_rot_init ≈ E_kin_rot_final 89 90 91 @info "Testing kinetic energy conservation with Two-term Taylor scheme" 92 sim = deepcopy(sim_init) 93 Granular.setTimeStep!(sim, epsilon=0.007) 94 tol = 0.02 95 @info "Relative tolerance: $(tol*100.)%" 96 Granular.run!(sim, temporal_integration_method="Two-term Taylor", verbose=verbose) 97 98 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim) 99 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim) 100 @test E_kin_lin_init ≈ E_kin_lin_final atol=E_kin_lin_init*tol 101 @test E_kin_rot_init ≈ E_kin_rot_final 102 103 104 @info "Testing kinetic energy conservation with Three-term Taylor scheme" 105 sim = deepcopy(sim_init) 106 Granular.setTimeStep!(sim, epsilon=0.07) 107 tol = 0.01 108 @info "Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)" 109 Granular.run!(sim, temporal_integration_method="Three-term Taylor", 110 verbose=verbose) 111 112 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim) 113 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim) 114 @test E_kin_lin_init ≈ E_kin_lin_final atol=E_kin_lin_init*tol 115 @test E_kin_rot_init ≈ E_kin_rot_final 116 117 118 @info "# Adding contact-normal viscosity" 119 @info "# One ice floe fixed" 120 sim = Granular.createSimulation(id="test") 121 Granular.addGrainCylindrical!(sim, [0., 0.], 10., 1., verbose=verbose) 122 Granular.addGrainCylindrical!(sim, [20.05, 0.], 10., 1., verbose=verbose) 123 sim.grains[1].lin_vel[1] = 0.1 124 sim.grains[1].contact_viscosity_normal = 1e4 125 sim.grains[2].contact_viscosity_normal = 1e4 126 sim.grains[2].fixed = true 127 128 E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim) 129 E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim) 130 131 # With decreasing timestep (epsilon towards 0), the explicit integration scheme 132 # should become more correct 133 134 Granular.setTotalTime!(sim, 10.0) 135 sim_init = deepcopy(sim) 136 137 138 @info "Testing kinetic energy conservation with Two-term Taylor scheme" 139 sim = deepcopy(sim_init) 140 Granular.setTimeStep!(sim, epsilon=0.007) 141 tol = 0.02 142 @info "Relative tolerance: $(tol*100.)%" 143 Granular.run!(sim, temporal_integration_method="Two-term Taylor", verbose=verbose) 144 145 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim) 146 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim) 147 @test E_kin_lin_init > E_kin_lin_final 148 @test E_kin_rot_init ≈ E_kin_rot_final 149 150 151 @info "Testing kinetic energy conservation with Three-term Taylor scheme" 152 sim = deepcopy(sim_init) 153 Granular.setTimeStep!(sim, epsilon=0.07) 154 tol = 0.01 155 @info "Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)" 156 Granular.run!(sim, temporal_integration_method="Three-term Taylor", 157 verbose=verbose) 158 159 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim) 160 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim) 161 @test E_kin_lin_init > E_kin_lin_final 162 @test E_kin_rot_init ≈ E_kin_rot_final 163 164 165 @info "# Ice floes free to move" 166 167 sim = Granular.createSimulation(id="test") 168 Granular.addGrainCylindrical!(sim, [0., 0.], 10., 1., verbose=verbose) 169 Granular.addGrainCylindrical!(sim, [20.05, 0.], 10., 1., verbose=verbose) 170 sim.grains[1].lin_vel[1] = 0.1 171 sim.grains[1].contact_viscosity_normal = 1e4 172 sim.grains[2].contact_viscosity_normal = 1e4 173 174 E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim) 175 E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim) 176 177 # With decreasing timestep (epsilon towards 0), the explicit integration scheme 178 # should become more correct 179 180 Granular.setTotalTime!(sim, 10.0) 181 sim_init = deepcopy(sim) 182 183 @info "Testing kinetic energy conservation with Two-term Taylor scheme" 184 sim = deepcopy(sim_init) 185 Granular.setTimeStep!(sim, epsilon=0.007) 186 tol = 0.02 187 @info "Relative tolerance: $(tol*100.)%" 188 Granular.run!(sim, temporal_integration_method="Two-term Taylor", verbose=verbose) 189 190 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim) 191 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim) 192 @test E_kin_lin_init > E_kin_lin_final 193 @test E_kin_rot_init ≈ E_kin_rot_final 194 195 196 @info "Testing kinetic energy conservation with Three-term Taylor scheme" 197 sim = deepcopy(sim_init) 198 Granular.setTimeStep!(sim, epsilon=0.07) 199 tol = 0.01 200 @info "Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)" 201 Granular.run!(sim, temporal_integration_method="Three-term Taylor", 202 verbose=verbose) 203 204 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim) 205 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim) 206 @test E_kin_lin_init > E_kin_lin_final 207 @test E_kin_rot_init ≈ E_kin_rot_final 208 209 210 @info "# Testing allow_*_acc for fixed grains" 211 sim = Granular.createSimulation(id="test") 212 Granular.addGrainCylindrical!(sim, [0., 0.], 10., 1., verbose=verbose) 213 Granular.addGrainCylindrical!(sim, [20.05, 0.], 10., 1., verbose=verbose) 214 sim.grains[1].lin_vel[1] = 0.1 215 sim.grains[2].fixed = true 216 217 E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim) 218 E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim) 219 grain2_pos_init = sim.grains[2].lin_pos 220 221 Granular.setTotalTime!(sim, 10.0) 222 Granular.setTimeStep!(sim, epsilon=0.07) 223 sim_init = deepcopy(sim) 224 sim.grains[2].allow_y_acc = true # should not influence result 225 226 @info "Two-term Taylor scheme: allow_y_acc" 227 sim = deepcopy(sim_init) 228 sim.grains[2].allow_y_acc = true # should not influence result 229 tol = 0.2 230 Granular.run!(sim, temporal_integration_method="Two-term Taylor", verbose=verbose) 231 232 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim) 233 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim) 234 @test E_kin_lin_init ≈ E_kin_lin_final atol=E_kin_lin_init*tol 235 @test E_kin_rot_init ≈ E_kin_rot_final 236 @test sim.grains[2].lin_pos ≈ grain2_pos_init 237 238 @info "Two-term Taylor scheme: allow_x_acc" 239 sim = deepcopy(sim_init) 240 sim.grains[2].allow_x_acc = true # should influence result 241 tol = 0.2 242 Granular.run!(sim, temporal_integration_method="Two-term Taylor", verbose=verbose) 243 244 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim) 245 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim) 246 @test E_kin_lin_init ≈ E_kin_lin_final atol=E_kin_lin_init*tol 247 @test E_kin_rot_init ≈ E_kin_rot_final 248 @test sim.grains[2].lin_pos[1] > grain2_pos_init[1] 249 250 @info "Three-term Taylor scheme: allow_y_acc" 251 sim = deepcopy(sim_init) 252 tol = 0.02 253 sim.grains[2].allow_y_acc = true # should influence result 254 Granular.run!(sim, temporal_integration_method="Three-term Taylor", verbose=verbose) 255 256 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim) 257 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim) 258 @test E_kin_lin_init ≈ E_kin_lin_final atol=E_kin_lin_init*tol 259 @test E_kin_rot_init ≈ E_kin_rot_final 260 @test sim.grains[2].lin_pos ≈ grain2_pos_init 261 262 @info "Three-term Taylor scheme: allow_x_acc" 263 sim = deepcopy(sim_init) 264 tol = 0.02 265 sim.grains[2].allow_x_acc = true # should influence result 266 Granular.run!(sim, temporal_integration_method="Three-term Taylor", verbose=verbose) 267 268 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim) 269 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim) 270 @test E_kin_lin_init ≈ E_kin_lin_final atol=E_kin_lin_init*tol 271 @test E_kin_rot_init ≈ E_kin_rot_final 272 @test sim.grains[2].lin_pos[1] > grain2_pos_init[1] 273 274 #= 275 @info "# Test stability under collision with fixed particles different allow_*_acc" 276 r = 10. 277 i = 1 278 for tensile_strength in [0.0, 200e3] 279 for angle in range(0, 2π, 7) 280 for allow_x_acc in [false, true] 281 for allow_y_acc in [false, true] 282 @info "Test $i" 283 @info "Contact angle: $angle rad" 284 @info "allow_x_acc = $allow_x_acc" 285 @info "allow_y_acc = $allow_y_acc" 286 @info "tensile_strength = $tensile_strength Pa" 287 288 sim = Granular.createSimulation() 289 sim.id = "test-$i-$allow_x_acc-$allow_y_acc-C=$tensile_strength" 290 Granular.addGrainCylindrical!(sim, [0., 0.], r, 1., verbose=verbose) 291 Granular.addGrainCylindrical!(sim, [2.0*r*cos(angle), 2.0*r*sin(angle)], 292 r, 1., verbose=verbose) 293 sim.grains[1].lin_vel = r/10.0 .* [cos(angle), sin(angle)] 294 295 E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim) 296 E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim) 297 grain1_pos_init = sim.grains[1].lin_pos 298 grain2_pos_init = sim.grains[2].lin_pos 299 300 sim.grains[1].fixed = true 301 sim.grains[2].fixed = true 302 303 sim.grains[1].allow_x_acc = allow_x_acc 304 sim.grains[2].allow_x_acc = allow_x_acc 305 sim.grains[1].allow_y_acc = allow_y_acc 306 sim.grains[2].allow_y_acc = allow_y_acc 307 308 sim.grains[1].tensile_strength = tensile_strength 309 sim.grains[2].tensile_strength = tensile_strength 310 311 Granular.setTotalTime!(sim, 20.0) 312 Granular.setTimeStep!(sim, epsilon=0.07) 313 sim_init = deepcopy(sim) 314 315 @info "TY3" 316 sim = deepcopy(sim_init) 317 tol = 0.02 318 Granular.setOutputFileInterval!(sim, 1.0) 319 Granular.run!(sim, temporal_integration_method="Three-term Taylor", 320 verbose=verbose) 321 Granular.render(sim) 322 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim) 323 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim) 324 @test E_kin_lin_init ≈ E_kin_lin_final atol=E_kin_lin_init*tol 325 326 i += 1 327 end 328 end 329 end 330 end 331 =#