tAdd additional bond tests and fix bond problems - 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
---
(DIR) commit 2ab3e5cc94c45a4baff6d1e5ca769c3eeebe0309
(DIR) parent c0fa76fc348cd9d1285191817827717f951e26cb
(HTM) Author: Anders Damsgaard <andersd@riseup.net>
Date: Thu, 22 Feb 2018 15:40:28 -0500
Add additional bond tests and fix bond problems
Diffstat:
M src/interaction.jl | 13 ++++++++-----
M test/cohesion.jl | 143 +++++++++++++++++++++++++++++--
2 files changed, 142 insertions(+), 14 deletions(-)
---
(DIR) diff --git a/src/interaction.jl b/src/interaction.jl
t@@ -238,12 +238,12 @@ function interactGrains!(simulation::Simulation, i::Int, j::Int, ic::Int)
if tensile_strength > 0.0
# Determine bending momentum on contact [N*m],
# (converting k_n to E to bar(k_n))
- M_t = -(k_n*R_ij/(A_ij*(simulation.grains[i].contact_radius +
- simulation.grains[j].contact_radius)))*I_ij*θ_t
+ M_t = (k_n*R_ij/(A_ij*(simulation.grains[i].contact_radius +
+ simulation.grains[j].contact_radius)))*I_ij*θ_t
end
# Reset contact age (breaking bond) if bond strength is exceeded
- if δ_n > 0.0 && abs(force_n)/A_ij + abs(M_t)*R_ij/I_ij > tensile_strength
+ if δ_n >= 0.0 && abs(force_n)/A_ij + abs(M_t)*R_ij/I_ij > tensile_strength
force_n = 0.
force_t = 0.
simulation.grains[i].contacts[ic] = 0 # remove contact
t@@ -252,7 +252,7 @@ function interactGrains!(simulation::Simulation, i::Int, j::Int, ic::Int)
end
# Limit compressive stress if the prefactor is set to a positive value
- if δ_n < 0.0 && compressive_strength > 0. &&
+ if δ_n <= 0.0 && compressive_strength > 0. &&
abs(force_n) >= compressive_strength
# Determine the overlap distance where yeild stress is reached
t@@ -311,7 +311,10 @@ function interactGrains!(simulation::Simulation, i::Int, j::Int, ic::Int)
error("unknown contact_tangential_rheology (k_t = $k_t, γ_t = $γ_t")
end
- if shear_strength > 0.0 && abs(force_t)/A > shear_strength
+ # Break bond under extension through bending failure
+ if δ_n < 0.0 && tensile_strength > 0.0 && shear_strength > 0.0 &&
+ abs(force_t)/A_ij > shear_strength
+
force_n = 0.
force_t = 0.
simulation.grains[i].contacts[ic] = 0 # remove contact
(DIR) diff --git a/test/cohesion.jl b/test/cohesion.jl
t@@ -38,17 +38,142 @@ Granular.run!(sim, verbose=verbose)
@test sim.grains[1].ang_vel ≈ 0.
@test sim.grains[2].ang_vel ≈ 0.
-info("# Add shear strength")
+info("# Add shear strength and test bending resistance (one grain rotating)")
sim = Granular.createSimulation(id="cohesion")
-Granular.addGrainCylindrical!(sim, [0., 0.], 10., 1., tensile_strength=500e3)
-Granular.addGrainCylindrical!(sim, [20.1, 0.], 10., 1., tensile_strength=500e3)
-sim.grains[1].lin_vel[1] = 0.1
+Granular.addGrainCylindrical!(sim, [0., 0.], 10.1, 1., tensile_strength=500e3,
+ shear_strength=500e3)
+Granular.addGrainCylindrical!(sim, [20., 0.], 10., 1., tensile_strength=500e3,
+ shear_strength=500e3)
+sim.grains[1].ang_vel = 0.1
+Granular.findContacts!(sim) # make sure contact is registered
+sim.grains[1].contact_radius=10.0 # decrease radius so there isn't compression
+E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim)
+E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim)
Granular.setTimeStep!(sim)
-Granular.setTotalTime!(sim, 10.)
+Granular.setTotalTime!(sim, 5.)
+#Granular.setOutputFileInterval!(sim, 0.1)
Granular.run!(sim, verbose=verbose)
-@test sim.grains[1].lin_vel[1] > 0.
+@test sim.grains[1].lin_vel[1] ≈ 0.
@test sim.grains[1].lin_vel[2] ≈ 0.
-@test sim.grains[2].lin_vel[1] > 0.
+@test sim.grains[2].lin_vel[1] ≈ 0.
@test sim.grains[2].lin_vel[2] ≈ 0.
-@test sim.grains[1].ang_vel ≈ 0.
-@test sim.grains[2].ang_vel ≈ 0.
+@test sim.grains[1].ang_vel != 0.
+@test sim.grains[2].ang_vel != 0.
+E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
+E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
+E_therm_final = Granular.totalGrainThermalEnergy(sim)
+@test E_kin_lin_init ≈ E_kin_lin_final
+@test E_kin_rot_init > E_kin_rot_final + E_therm_final
+
+info("# Add shear strength and test bending resistance (one grain rotating)")
+sim = Granular.createSimulation(id="cohesion")
+Granular.addGrainCylindrical!(sim, [0., 0.], 10.1, 1., tensile_strength=500e3,
+ shear_strength=500e3)
+Granular.addGrainCylindrical!(sim, [20., 0.], 10., 1., tensile_strength=500e3,
+ shear_strength=500e3)
+sim.grains[2].ang_vel = 0.1
+Granular.findContacts!(sim) # make sure contact is registered
+sim.grains[1].contact_radius=10.0 # decrease radius so there isn't compression
+E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim)
+E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim)
+Granular.setTimeStep!(sim)
+Granular.setTotalTime!(sim, 5.)
+#Granular.setOutputFileInterval!(sim, 0.1)
+Granular.run!(sim, verbose=verbose)
+@test sim.grains[1].lin_vel[1] ≈ 0.
+@test sim.grains[1].lin_vel[2] ≈ 0.
+@test sim.grains[2].lin_vel[1] ≈ 0.
+@test sim.grains[2].lin_vel[2] ≈ 0.
+@test sim.grains[1].ang_vel != 0.
+@test sim.grains[2].ang_vel != 0.
+E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
+E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
+E_therm_final = Granular.totalGrainThermalEnergy(sim)
+@test E_kin_lin_init ≈ E_kin_lin_final
+@test E_kin_rot_init > E_kin_rot_final + E_therm_final
+
+info("# Add shear strength and test bending resistance (both grains rotating)")
+sim = Granular.createSimulation(id="cohesion")
+Granular.addGrainCylindrical!(sim, [0., 0.], 10.0000001, 1., tensile_strength=500e3,
+ shear_strength=500e3)
+Granular.addGrainCylindrical!(sim, [20., 0.], 10., 1., tensile_strength=500e3,
+ shear_strength=500e3)
+sim.grains[1].ang_vel = 0.1
+sim.grains[2].ang_vel = -0.1
+Granular.findContacts!(sim) # make sure contact is registered
+sim.grains[1].contact_radius=10.0 # decrease radius so there isn't compression
+E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim)
+E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim)
+Granular.setTimeStep!(sim)
+Granular.setTotalTime!(sim, 5.)
+#Granular.setOutputFileInterval!(sim, 0.1)
+Granular.run!(sim, verbose=verbose)
+@test sim.grains[1].lin_vel[1] ≈ 0.
+@test sim.grains[1].lin_vel[2] ≈ 0.
+@test sim.grains[2].lin_vel[1] ≈ 0.
+@test sim.grains[2].lin_vel[2] ≈ 0.
+@test sim.grains[1].ang_vel != 0.
+@test sim.grains[2].ang_vel != 0.
+E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
+E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
+E_therm_final = Granular.totalGrainThermalEnergy(sim)
+@test E_kin_lin_init ≈ E_kin_lin_final
+@test E_kin_rot_init > E_kin_rot_final + E_therm_final
+
+info("# Break bond through bending I")
+sim = Granular.createSimulation(id="cohesion")
+Granular.addGrainCylindrical!(sim, [0., 0.], 10.0000001, 1., tensile_strength=500e3,
+ shear_strength=500e3)
+Granular.addGrainCylindrical!(sim, [20., 0.], 10., 1., tensile_strength=500e3,
+ shear_strength=500e3)
+sim.grains[1].ang_vel = 100
+sim.grains[2].ang_vel = -100
+Granular.findContacts!(sim) # make sure contact is registered
+sim.grains[1].contact_radius=10.0 # decrease radius so there isn't compression
+E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim)
+E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim)
+Granular.setTimeStep!(sim)
+Granular.setTotalTime!(sim, 5.)
+#Granular.setOutputFileInterval!(sim, 0.1)
+Granular.run!(sim, verbose=verbose)
+@test sim.grains[1].lin_vel[1] ≈ 0.
+@test sim.grains[1].lin_vel[2] ≈ 0.
+@test sim.grains[2].lin_vel[1] ≈ 0.
+@test sim.grains[2].lin_vel[2] ≈ 0.
+@test sim.grains[1].ang_vel != 0.
+@test sim.grains[2].ang_vel != 0.
+E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
+E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
+E_therm_final = Granular.totalGrainThermalEnergy(sim)
+@test E_kin_lin_init ≈ E_kin_lin_final
+@test sim.grains[1].n_contacts == 0
+@test sim.grains[2].n_contacts == 0
+
+info("# Break bond through bending II")
+sim = Granular.createSimulation(id="cohesion")
+Granular.addGrainCylindrical!(sim, [0., 0.], 10.1, 1., tensile_strength=500e3,
+ shear_strength=50e3)
+Granular.addGrainCylindrical!(sim, [20., 0.], 10., 1., tensile_strength=500e3,
+ shear_strength=50e3)
+sim.grains[1].ang_vel = 100
+sim.grains[2].ang_vel = 0.0
+Granular.findContacts!(sim) # make sure contact is registered
+sim.grains[1].contact_radius=10.0 # decrease radius so there isn't compression
+E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim)
+E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim)
+Granular.setTimeStep!(sim)
+Granular.setTotalTime!(sim, 5.)
+#Granular.setOutputFileInterval!(sim, 0.1)
+Granular.run!(sim, verbose=verbose)
+@test sim.grains[1].lin_vel[1] ≈ 0.
+@test sim.grains[1].lin_vel[2] ≈ 0.
+@test sim.grains[2].lin_vel[1] ≈ 0.
+@test sim.grains[2].lin_vel[2] ≈ 0.
+@test sim.grains[1].ang_vel != 0.
+@test sim.grains[2].ang_vel != 0.
+E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
+E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
+E_therm_final = Granular.totalGrainThermalEnergy(sim)
+@test E_kin_lin_init ≈ E_kin_lin_final
+@test sim.grains[1].n_contacts == 0
+@test sim.grains[2].n_contacts == 0