cohesion.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