tAdd untested parameterization including bond bending and shear strength - 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 c0fa76fc348cd9d1285191817827717f951e26cb
(DIR) parent ee17a2d649961009a9f878201878098572f64596
(HTM) Author: Anders Damsgaard <andersd@riseup.net>
Date: Thu, 22 Feb 2018 13:43:35 -0500
Add untested parameterization including bond bending and shear strength
Diffstat:
M src/contact_search.jl | 1 +
M src/datatypes.jl | 7 +++++--
M src/grain.jl | 40 ++++++++++++++++++++++----------
M src/interaction.jl | 119 ++++++++++++++++++-------------
M src/io.jl | 32 ++++++++++++++++++++++++-------
M test/cohesion.jl | 15 +++++++++++++++
M test/vtk.jl | 4 ++--
7 files changed, 147 insertions(+), 71 deletions(-)
---
(DIR) diff --git a/src/contact_search.jl b/src/contact_search.jl
t@@ -323,6 +323,7 @@ function checkAndAddContact!(sim::Simulation, i::Int, j::Int,
position_ij
@inbounds fill!(sim.grains[i].
contact_parallel_displacement[ic] , 0.)
+ @inbounds sim.grains[i].contact_rotation[ic] = 0.
@inbounds sim.grains[i].contact_age[ic] = 0.
break
end
(DIR) diff --git a/src/datatypes.jl b/src/datatypes.jl
t@@ -47,7 +47,8 @@ mutable struct GrainCylindrical
youngs_modulus::Float64
poissons_ratio::Float64
tensile_strength::Float64
- tensile_heal_rate::Float64
+ shear_strength::Float64
+ strength_heal_rate::Float64
fracture_toughness::Float64
# Ocean/atmosphere interaction parameters
t@@ -64,6 +65,7 @@ mutable struct GrainCylindrical
contacts::Vector{Int}
position_vector::Vector{Vector{Float64}}
contact_parallel_displacement::Vector{Vector{Float64}}
+ contact_rotation::Vector{Float64}
contact_age::Vector{Float64}
granular_stress::Vector{Float64}
t@@ -140,7 +142,8 @@ mutable struct GrainArrays
youngs_modulus::Vector{Float64}
poissons_ratio::Vector{Float64}
tensile_strength::Vector{Float64}
- tensile_heal_rate::Vector{Float64}
+ shear_strength::Vector{Float64}
+ strength_heal_rate::Vector{Float64}
fracture_toughness::Vector{Float64}
ocean_drag_coeff_vert::Vector{Float64}
(DIR) diff --git a/src/grain.jl b/src/grain.jl
t@@ -14,7 +14,8 @@ export addGrainCylindrical!
contact_static_friction,
contact_dynamic_friction,
youngs_modulus, poissons_ratio,
- tensile_strength, tensile_heal_rate,
+ tensile_strength, shear_strength,
+ strength_heal_rate,
fracture_toughness,
ocean_drag_coeff_vert,
ocean_drag_coeff_horiz,
t@@ -69,8 +70,9 @@ are optional, and come with default values. The only required arguments are
contact-tangential stiffness from `youngs_modulus` [-].
* `tensile_strength::Float64 = 0.`: contact-tensile (cohesive) bond strength
[Pa].
-* `tensile_heal_rate::Float64 = 0.`: rate at which contact-tensile bond strength
- is obtained [Pa/s].
+* `shear_strength::Float64 = 0.`: shear strength of bonded contacts [Pa].
+* `strength_heal_rate::Float64 = 0.`: rate at which contact bond
+ strength is obtained [Pa/s].
* `fracture_toughness::Float64 = 0.`: maximum compressive
strength on granular contact (not currently enforced) [m^{1/2}*Pa]. A value
of 1.285e3 m^{1/2}*Pa is used for sea ice by Hopkins 2004.
t@@ -114,12 +116,13 @@ meter:
Granular.addGrainCylindrical!(sim, [1., 2.], 1., .5)
```
-The following example will create a grain with tensile strength (cohesion),
-and a velocity of 0.5 m/s towards -x:
+The following example will create a grain with tensile and shear strength, and a
+velocity of 0.5 m/s towards -x:
```julia
Granular.addGrainCylindrical!(sim, [4., 2.], 1., .5,
tensile_strength = 200e3,
+ shear_strength = 100e3,
lin_vel = [-.5, 0.])
```
t@@ -154,7 +157,8 @@ function addGrainCylindrical!(simulation::Simulation,
youngs_modulus::Float64 = 2e7,
poissons_ratio::Float64 = 0.185, # Hopkins 2004
tensile_strength::Float64 = 0.,
- tensile_heal_rate::Float64 = Inf,
+ shear_strength::Float64 = 0.,
+ strength_heal_rate::Float64 = Inf,
fracture_toughness::Float64 = 0.,
ocean_drag_coeff_vert::Float64 = 0.85, # H&C 2011
ocean_drag_coeff_horiz::Float64 = 5e-4, # H&C 2011
t@@ -203,6 +207,7 @@ function addGrainCylindrical!(simulation::Simulation,
contacts::Array{Int, 1} = zeros(Int, simulation.Nc_max)
position_vector = Vector{Vector{Float64}}(simulation.Nc_max)
contact_parallel_displacement = Vector{Vector{Float64}}(simulation.Nc_max)
+ contact_rotation::Vector{Float64} = zeros(Float64, simulation.Nc_max)
contact_age::Vector{Float64} = zeros(Float64, simulation.Nc_max)
for i=1:simulation.Nc_max
position_vector[i] = zeros(2)
t@@ -250,7 +255,8 @@ function addGrainCylindrical!(simulation::Simulation,
youngs_modulus,
poissons_ratio,
tensile_strength,
- tensile_heal_rate,
+ shear_strength,
+ strength_heal_rate,
fracture_toughness,
ocean_drag_coeff_vert,
t@@ -265,6 +271,7 @@ function addGrainCylindrical!(simulation::Simulation,
contacts,
position_vector,
contact_parallel_displacement,
+ contact_rotation,
contact_age,
granular_stress,
t@@ -413,7 +420,9 @@ function convertGrainDataToArrays(simulation::Simulation)
Array{Float64}(length(simulation.grains)),
## tensile_strength
Array{Float64}(length(simulation.grains)),
- ## tensile_heal_rate
+ ## shear_strength
+ Array{Float64}(length(simulation.grains)),
+ ## strength_heal_rate
Array{Float64}(length(simulation.grains)),
## fracture_toughness
Array{Float64}(length(simulation.grains)),
t@@ -500,7 +509,8 @@ function convertGrainDataToArrays(simulation::Simulation)
ifarr.youngs_modulus[i] = simulation.grains[i].youngs_modulus
ifarr.poissons_ratio[i] = simulation.grains[i].poissons_ratio
ifarr.tensile_strength[i] = simulation.grains[i].tensile_strength
- ifarr.tensile_heal_rate[i] = simulation.grains[i].tensile_heal_rate
+ ifarr.shear_strength[i] = simulation.grains[i].shear_strength
+ ifarr.strength_heal_rate[i] = simulation.grains[i].strength_heal_rate
ifarr.fracture_toughness[i] =
simulation.grains[i].fracture_toughness
t@@ -574,7 +584,8 @@ function deleteGrainArrays!(ifarr::GrainArrays)
ifarr.youngs_modulus = f1
ifarr.poissons_ratio = f1
ifarr.tensile_strength = f1
- ifarr.tensile_heal_rate = f1
+ ifarr.shear_strength = f1
+ ifarr.strength_heal_rate = f1
ifarr.fracture_toughness = f1
ifarr.ocean_drag_coeff_vert = f1
t@@ -641,6 +652,8 @@ function printGrainInfo(f::GrainCylindrical)
println(" E: $(f.youngs_modulus) Pa")
println(" ν: $(f.poissons_ratio)")
println(" tensile_strength: $(f.tensile_strength) Pa")
+ println(" shear_strength: $(f.shear_strength) Pa")
+ println(" strength_heal_rate: $(f.strength_heal_rate) Pa/s")
println(" fracture_toughness: $(f.fracture_toughness) m^0.5 Pa\n")
println(" c_o_v: $(f.ocean_drag_coeff_vert)")
t@@ -652,7 +665,8 @@ function printGrainInfo(f::GrainCylindrical)
println(" n_contacts: $(f.n_contacts)")
println(" contacts: $(f.contacts)")
println(" pos_ij: $(f.position_vector)\n")
- println(" delta_t: $(f.contact_parallel_displacement)\n")
+ println(" δ_t: $(f.contact_parallel_displacement)\n")
+ println(" θ_t: $(f.contact_rotation)\n")
println(" granular_stress: $(f.granular_stress) Pa")
println(" ocean_stress: $(f.ocean_stress) Pa")
t@@ -802,7 +816,8 @@ function compareGrains(if1::GrainCylindrical, if2::GrainCylindrical)
@test if1.youngs_modulus ≈ if2.youngs_modulus
@test if1.poissons_ratio ≈ if2.poissons_ratio
@test if1.tensile_strength ≈ if2.tensile_strength
- @test if1.tensile_heal_rate ≈ if2.tensile_heal_rate
+ @test if1.shear_strength ≈ if2.shear_strength
+ @test if1.strength_heal_rate ≈ if2.strength_heal_rate
@test if1.fracture_toughness ≈
if2.fracture_toughness
t@@ -820,6 +835,7 @@ function compareGrains(if1::GrainCylindrical, if2::GrainCylindrical)
@test if1.position_vector == if2.position_vector
@test if1.contact_parallel_displacement ==
if2.contact_parallel_displacement
+ @test if1.contact_rotation == if2.contact_rotation
@test if1.contact_age ≈ if2.contact_age
@test if1.granular_stress ≈ if2.granular_stress
(DIR) diff --git a/src/interaction.jl b/src/interaction.jl
t@@ -125,46 +125,51 @@ function interactGrains!(simulation::Simulation, i::Int, j::Int, ic::Int)
# Inter-position vector, use stored value which is corrected for boundary
# periodicity
- p = simulation.grains[i].position_vector[ic]
- dist = norm(p)
+ const p = simulation.grains[i].position_vector[ic]
+ const dist = norm(p)
- r_i = simulation.grains[i].contact_radius
- r_j = simulation.grains[j].contact_radius
+ const r_i = simulation.grains[i].contact_radius
+ const r_j = simulation.grains[j].contact_radius
# Floe distance; <0: compression, >0: tension
- δ_n = dist - (r_i + r_j)
+ const δ_n = dist - (r_i + r_j)
# Local axes
- n = p/dist
- t = [-n[2], n[1]]
+ const n = p/dist
+ const t = [-n[2], n[1]]
# Contact kinematics
- vel_lin = simulation.grains[i].lin_vel -
+ const vel_lin = simulation.grains[i].lin_vel -
simulation.grains[j].lin_vel
- vel_n = dot(vel_lin, n)
- vel_t = dot(vel_lin, t) -
+ const vel_n = dot(vel_lin, n)
+ const vel_t = dot(vel_lin, t) -
harmonicMean(r_i, r_j)*(simulation.grains[i].ang_vel +
simulation.grains[j].ang_vel)
# Correct old tangential displacement for contact rotation and add new
- δ_t0 = simulation.grains[i].contact_parallel_displacement[ic]
- δ_t = dot(t, δ_t0 - (n*dot(n, δ_t0))) + vel_t*simulation.time_step
+ const δ_t0 = simulation.grains[i].contact_parallel_displacement[ic]
+ const δ_t = dot(t, δ_t0 - (n*dot(n, δ_t0))) + vel_t*simulation.time_step
+
+ # Determine the contact rotation
+ const θ_t = simulation.grains[i].contact_rotation[ic] +
+ (simulation.grains[j].ang_vel - simulation.grains[i].ang_vel) *
+ simulation.time_step
# Effective radius
- R_ij = harmonicMean(r_i, r_j)
+ const R_ij = harmonicMean(r_i, r_j)
# Contact area
- A_ij = R_ij*min(simulation.grains[i].thickness,
- simulation.grains[j].thickness)
+ const A_ij = R_ij*min(simulation.grains[i].thickness,
+ simulation.grains[j].thickness)
# Contact mechanical parameters
if simulation.grains[i].youngs_modulus > 0. &&
simulation.grains[j].youngs_modulus > 0.
- E = harmonicMean(simulation.grains[i].youngs_modulus,
- simulation.grains[j].youngs_modulus)
- ν = harmonicMean(simulation.grains[i].poissons_ratio,
- simulation.grains[j].poissons_ratio)
+ const E = harmonicMean(simulation.grains[i].youngs_modulus,
+ simulation.grains[j].youngs_modulus)
+ const ν = harmonicMean(simulation.grains[i].poissons_ratio,
+ simulation.grains[j].poissons_ratio)
# Effective normal and tangential stiffness
k_n = E * A_ij/R_ij
t@@ -179,14 +184,14 @@ function interactGrains!(simulation::Simulation, i::Int, j::Int, ic::Int)
simulation.grains[j].contact_stiffness_tangential)
end
- γ_n = harmonicMean(simulation.grains[i].contact_viscosity_normal,
- simulation.grains[j].contact_viscosity_normal)
+ const γ_n = harmonicMean(simulation.grains[i].contact_viscosity_normal,
+ simulation.grains[j].contact_viscosity_normal)
- γ_t = harmonicMean(simulation.grains[i].contact_viscosity_tangential,
- simulation.grains[j].contact_viscosity_tangential)
+ const γ_t = harmonicMean(simulation.grains[i].contact_viscosity_tangential,
+ simulation.grains[j].contact_viscosity_tangential)
- μ_d_minimum = min(simulation.grains[i].contact_dynamic_friction,
- simulation.grains[j].contact_dynamic_friction)
+ const μ_d_minimum = min(simulation.grains[i].contact_dynamic_friction,
+ simulation.grains[j].contact_dynamic_friction)
# Determine contact forces
if k_n ≈ 0. && γ_n ≈ 0. # No interaction
t@@ -216,28 +221,39 @@ function interactGrains!(simulation::Simulation, i::Int, j::Int, ic::Int)
idx_weakest = j
end
- # Add tensile strength during extension or limit compressive strength
- if δ_n > 0.
- # Contact tensile strength increases linearly with contact age until
- # tensile stress exceeds tensile strength
-
- # linearly increase tensile strength with time until max. value
- tensile_strength = min(simulation.grains[i].contact_age[ic]*
- simulation.grains[i].tensile_heal_rate,
- simulation.grains[i].tensile_strength)
-
+ # Grain-pair moment of inertia [m^4]
+ const I_ij = 2.0/3.0*R_ij^3*min(simulation.grains[i].thickness,
+ simulation.grains[j].thickness)
+
+
+ # Contact tensile strength increases linearly with contact age until
+ # tensile stress exceeds tensile strength.
+ tensile_strength = min(simulation.grains[i].contact_age[ic]*
+ simulation.grains[i].strength_heal_rate,
+ simulation.grains[i].tensile_strength)
+ shear_strength = min(simulation.grains[i].contact_age[ic]*
+ simulation.grains[i].strength_heal_rate,
+ simulation.grains[i].shear_strength)
+ M_t = 0.0
+ 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
+ end
- # break bond
- if abs(force_n) >= tensile_strength*A_ij
- force_n = 0.
- force_t = 0.
- simulation.grains[i].contacts[ic] = 0 # remove contact
- simulation.grains[i].n_contacts -= 1
- simulation.grains[j].n_contacts -= 1
- 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
+ force_n = 0.
+ force_t = 0.
+ simulation.grains[i].contacts[ic] = 0 # remove contact
+ simulation.grains[i].n_contacts -= 1
+ simulation.grains[j].n_contacts -= 1
+ end
# Limit compressive stress if the prefactor is set to a positive value
- elseif compressive_strength > 0. && abs(force_n) >= compressive_strength
+ if δ_n < 0.0 && compressive_strength > 0. &&
+ abs(force_n) >= compressive_strength
# Determine the overlap distance where yeild stress is reached
δ_n_yield = -compressive_strength*A_ij/k_n
t@@ -260,7 +276,6 @@ function interactGrains!(simulation::Simulation, i::Int, j::Int, ic::Int)
# Coulomb slip
if force_t > μ_d_minimum*abs(force_n)
force_t = μ_d_minimum*abs(force_n)
- simulation.grains[i].contact_age[ic] = 0.0
# Nguyen et al 2009 "Thermomechanical modelling of friction effects
# in granular flows using the discrete element method"
t@@ -282,7 +297,6 @@ function interactGrains!(simulation::Simulation, i::Int, j::Int, ic::Int)
if abs(force_t) > μ_d_minimum*abs(force_n)
force_t = μ_d_minimum*abs(force_n)*force_t/abs(force_t)
δ_t = (-force_t - γ_t*vel_t)/k_t
- simulation.grains[i].contact_age[ic] = 0.0
# Nguyen et al 2009 "Thermomechanical modelling of friction effects
# in granular flows using the discrete element method"
t@@ -297,14 +311,23 @@ 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
+ force_n = 0.
+ force_t = 0.
+ simulation.grains[i].contacts[ic] = 0 # remove contact
+ simulation.grains[i].n_contacts -= 1
+ simulation.grains[j].n_contacts -= 1
+ end
+
simulation.grains[i].contact_parallel_displacement[ic] = δ_t*t
+ simulation.grains[i].contact_rotation[ic] = θ_t
simulation.grains[i].contact_age[ic] += simulation.time_step
simulation.grains[i].force += force_n*n + force_t*t;
simulation.grains[j].force -= force_n*n + force_t*t;
- simulation.grains[i].torque += -force_t*R_ij
- simulation.grains[j].torque += -force_t*R_ij
+ simulation.grains[i].torque += -force_t*R_ij + M_t
+ simulation.grains[j].torque += -force_t*R_ij - M_t
simulation.grains[i].pressure +=
force_n/simulation.grains[i].side_surface_area;
(DIR) diff --git a/src/io.jl b/src/io.jl
t@@ -402,8 +402,10 @@ function writeGrainVTK(simulation::Simulation,
"Poisson's ratio [-]")
WriteVTK.vtk_point_data(vtkfile, ifarr.tensile_strength,
"Tensile strength [Pa]")
- WriteVTK.vtk_point_data(vtkfile, ifarr.tensile_heal_rate,
- "Tensile healing rate [1/s]")
+ WriteVTK.vtk_point_data(vtkfile, ifarr.shear_strength,
+ "Shear strength [Pa]")
+ WriteVTK.vtk_point_data(vtkfile, ifarr.strength_heal_rate,
+ "Bond strength healing rate [Pa/s]")
WriteVTK.vtk_point_data(vtkfile, ifarr.fracture_toughness,
"Fracture toughness [m^0.5 Pa]")
t@@ -467,6 +469,7 @@ function writeGrainInteractionVTK(simulation::Simulation,
contact_stiffness = Float64[]
tensile_stress = Float64[]
shear_displacement = Vector{Float64}[]
+ contact_rotation = Float64[]
contact_age = Float64[]
for i=1:length(simulation.grains)
t@@ -528,8 +531,10 @@ function writeGrainInteractionVTK(simulation::Simulation,
push!(contact_stiffness, k_n)
push!(tensile_stress, k_n*δ_n/A_ij)
- push!(shear_displacement, simulation.grains[i].
- contact_parallel_displacement[ic])
+ push!(shear_displacement,
+ simulation.grains[i]. contact_parallel_displacement[ic])
+ push!(contact_rotation,
+ simulation.grains[i].contact_rotation[ic])
push!(contact_age, simulation.grains[i].contact_age[ic])
end
t@@ -620,6 +625,15 @@ function writeGrainInteractionVTK(simulation::Simulation,
write(f, " </DataArray>\n")
write(f, " <DataArray type=\"Float32\" " *
+ "Name=\"Contact rotation [rad]\" NumberOfComponents=\"1\"
+ format=\"ascii\">\n")
+ for i=1:length(i1)
+ @inbounds write(f, "$(contact_rotation[i]) ")
+ end
+ write(f, "\n")
+ write(f, " </DataArray>\n")
+
+ write(f, " <DataArray type=\"Float32\" " *
"Name=\"Contact age [s]\" NumberOfComponents=\"1\"
format=\"ascii\">\n")
for i=1:length(i1)
t@@ -853,7 +867,8 @@ imagegrains.PointArrayStatus = [
"Young's modulus [Pa]",
"Poisson's ratio [-]",
'Tensile strength [Pa]'
-'Tensile heal rate [1/s]'
+'Shear strength [Pa]'
+'Strength heal rate [Pa/s]'
'Compressive strength prefactor [m^0.5 Pa]',
'Ocean drag coefficient (vertical) [-]',
'Ocean drag coefficient (horizontal) [-]',
t@@ -1278,6 +1293,7 @@ function plotGrains(sim::Simulation;
contact_stiffness = Float64[]
tensile_stress = Float64[]
shear_displacement = Vector{Float64}[]
+ contact_rotation = Float64[]
contact_age = Float64[]
for i=1:length(sim.grains)
for ic=1:sim.Nc_max
t@@ -1324,8 +1340,10 @@ function plotGrains(sim::Simulation;
push!(contact_stiffness, k_n)
push!(tensile_stress, k_n*δ_n/A_ij)
- push!(shear_displacement, sim.grains[i].
- contact_parallel_displacement[ic])
+ push!(shear_displacement,
+ sim.grains[i].contact_parallel_displacement[ic])
+ push!(contact_rotation,
+ sim.grains[i].contact_rotation[ic])
push!(contact_age, sim.grains[i].contact_age[ic])
end
(DIR) diff --git a/test/cohesion.jl b/test/cohesion.jl
t@@ -37,3 +37,18 @@ Granular.run!(sim, verbose=verbose)
@test sim.grains[2].lin_vel[2] ≈ 0.
@test sim.grains[1].ang_vel ≈ 0.
@test sim.grains[2].ang_vel ≈ 0.
+
+info("# Add shear strength")
+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.setTimeStep!(sim)
+Granular.setTotalTime!(sim, 10.)
+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.
(DIR) diff --git a/test/vtk.jl b/test/vtk.jl
t@@ -29,12 +29,12 @@ end
grainpath = "test/test.grains.1.vtu"
grainchecksum =
-"f8aa4ab923b2466d7a6ffb835bfff4d0e0ff0a168ff267431a082bb3018caaa1 " *
+"eff130f14975dd5da06186c5b61c0c1a9d679f2188d03019c8144898ddc902b6 " *
grainpath * "\n"
graininteractionpath = "test/test.grain-interaction.1.vtp"
graininteractionchecksum =
-"881598f8f7279ece4301f6c94cb8f9146eb695f8c710edb446f49c1f7a061b84 " *
+"c61f314a997c4405eaa98e6a4e3f81368ab2e1c60d06d9a0d998b69d7c8fb1bf " *
graininteractionpath * "\n"
oceanpath = "test/test.ocean.1.vts"