commit 5e16f290d31e8d0d24032fc60ec40a5dfbcf52f6
parent b3a11d3fab0cc139a56290b080d1b68917566ee8
Author: Anders Damsgaard <andersd@riseup.net>
Date: Tue, 9 May 2017 12:12:52 -0400
do not store overlap, begin implementing tangential elasticity
Diffstat:
5 files changed, 35 insertions(+), 40 deletions(-)
diff --git a/src/contact_search.jl b/src/contact_search.jl
@@ -136,7 +136,6 @@ function checkAndAddContact!(simulation::Simulation, i::Int, j::Int)
# Check if grains overlap (overlap when negative)
if overlap_ij < 0.0
push!(simulation.contact_pairs, [i, j])
- push!(simulation.overlaps, overlap_ij*position_ij/norm(position_ij))
push!(simulation.contact_parallel_displacement, zeros(2))
end
end
diff --git a/src/datatypes.jl b/src/datatypes.jl
@@ -172,7 +172,6 @@ type Simulation
ice_floes::Array{IceFloeCylindrical, 1}
contact_pairs::Array{Array{Int, 1}, 1}
- overlaps::Array{Array{Float64, 1}, 1}
contact_parallel_displacement::Array{Array{Float64, 1}, 1}
ocean::Ocean
diff --git a/src/interaction.jl b/src/interaction.jl
@@ -11,11 +11,10 @@ function interact!(simulation::Simulation;
# IceFloe to grain collisions
while !isempty(simulation.contact_pairs)
contact_pair = pop!(simulation.contact_pairs)
- overlap = pop!(simulation.overlaps)
contact_parallel_displacement =
pop!(simulation.contact_parallel_displacement)
interactIceFloes!(simulation, contact_pair[1], contact_pair[2],
- overlap, contact_parallel_displacement,
+ contact_parallel_displacement,
contact_normal_rheology=contact_normal_rheology,
contact_tangential_rheology=
contact_tangential_rheology)
@@ -30,23 +29,24 @@ function adds the compressive force of the interaction to the ice floe
"""
function interactIceFloes!(simulation::Simulation,
i::Int, j::Int,
- overlap::vector,
contact_parallel_displacement::vector;
contact_normal_rheology::String = "Linear Elastic",
contact_tangential_rheology::String = "None")
- force_n = 0.
- force_t = 0.
+ force_n = 0. # Contact-normal force
+ force_t = 0. # Contact-parallel (tangential) force
+ # Inter-position vector
p = simulation.ice_floes[i].lin_pos - simulation.ice_floes[j].lin_pos
dist = norm(p)
r_i = simulation.ice_floes[i].contact_radius
r_j = simulation.ice_floes[j].contact_radius
- dn = dist - (r_i + r_j)
+ # Floe distance
+ delta_n = dist - (r_i + r_j)
- if dn < 0. # Contact (this should always occur)
+ if delta_n < 0. # Contact (this should always occur)
# Local axes
n = p/dist
@@ -60,9 +60,15 @@ function interactIceFloes!(simulation::Simulation,
harmonicMean(r_i, r_j)*(simulation.ice_floes[i].ang_vel +
simulation.ice_floes[j].ang_vel)
+ # Correct old tangential displacement for contact rotation and add new
+ delta_t = dot(t, contact_parallel_displacement -
+ (n*dot(n, contact_parallel_displacement))) +
+ vel_t*simulation.time_step
+
+ # Effective radius
R_ij = harmonicMean(simulation.ice_floes[i].contact_radius,
- simulation.ice_floes[j].contact_radius) -
- norm(overlap)/2.
+ simulation.ice_floes[j].contact_radius) -
+ abs(delta_n)/2.
# Contact mechanical parameters
k_n_harmonic_mean =
@@ -87,29 +93,41 @@ function interactIceFloes!(simulation::Simulation,
# Determine contact forces
if contact_normal_rheology == "None"
force_n = 0.
+
elseif contact_normal_rheology == "Linear Elastic"
- force_n = -k_n_harmonic_mean*dn
+ force_n = -k_n_harmonic_mean*delta_n
+
elseif contact_normal_rheology == "Linear Viscous Elastic"
- force_n = -k_n_harmonic_mean*dn - gamma_n_harmonic_mean*vel_n
+ force_n = -k_n_harmonic_mean*delta_n - gamma_n_harmonic_mean*vel_n
if force_n < 0.
force_n = 0.
end
+
else
error("unknown contact_normal_rheology '$contact_normal_rheology'")
end
if contact_tangential_rheology == "None"
# do nothing
- elseif contact_tangential_rheology == "Linear Elastic Frictional"
- error("not yet implemented")
+
+ elseif contact_tangential_rheology ==
+ "Linear Elastic Viscous Frictional"
+ force_t = -k_t_harmonic_mean*delta_t - gamma_t_harmonic_mean*vel_t
+ if abs(force_t) > mu_d_minimum*abs(force_n)
+ force_t = mu_d_minimum*abs(force_n)*force_t/abs(force_t)
+ delta_t = (-force_t - gamma_t_harmonic_mean*vel_t)/
+ k_t_harmonic_mean
+ end
+
elseif contact_tangential_rheology == "Linear Viscous Frictional"
force_t = abs(gamma_t_harmonic_mean * vel_t)
- if force_t > mu_d_minimum*norm(force_n)
- force_t = mu_d_minimum*norm(force_n)
+ if force_t > mu_d_minimum*abs(force_n)
+ force_t = mu_d_minimum*abs(force_n)
end
if vel_t > 0.
force_t = -force_t
end
+
else
error("unknown contact_tangential_rheology ",
"'$contact_tangential_rheology'")
@@ -120,6 +138,8 @@ function interactIceFloes!(simulation::Simulation,
"floes $i and $j")
end
+ contact_parallel_displacement = delta_t*t
+
simulation.ice_floes[i].force += force_n*n + force_t*t;
simulation.ice_floes[j].force -= force_n*n + force_t*t;
diff --git a/src/simulation.jl b/src/simulation.jl
@@ -11,7 +11,6 @@ export createSimulation
file_number::Int=0,
ice_floes=Array{IceFloeCylindrical, 1}[],
contact_pairs=Array{Int64, 1}[],
- overlaps=Array{Array{Float64, 1}, 1}[],
contact_parallel_displacement=
Array{Array{Float64, 1}, 1}[])
@@ -31,7 +30,6 @@ function createSimulation(;id::String="unnamed",
file_time_since_output_file::Float64=0.,
ice_floes=Array{IceFloeCylindrical, 1}[],
contact_pairs=Array{Int64, 1}[],
- overlaps=Array{Array{Float64, 1}, 1}[],
contact_parallel_displacement=
Array{Array{Float64, 1}, 1}[],
ocean::Ocean=createEmptyOcean())
@@ -46,7 +44,6 @@ function createSimulation(;id::String="unnamed",
file_time_since_output_file,
ice_floes,
contact_pairs,
- overlaps,
contact_parallel_displacement,
ocean)
end
diff --git a/test/contact-search-and-geometry.jl b/test/contact-search-and-geometry.jl
@@ -20,10 +20,8 @@ info("Testing findContactsAllToAll(...)")
sim_copy = deepcopy(sim)
SeaIce.findContactsAllToAll!(sim)
-@test 1 == length(sim.overlaps)
@test 1 == length(sim.contact_pairs)
@test_approx_eq [1, 2] sim.contact_pairs[1]
-@test_approx_eq [2., 0.] sim.overlaps[1]
info("Testing findContacts(...)")
@@ -31,19 +29,15 @@ sim = deepcopy(sim_copy)
SeaIce.findContacts!(sim)
sim.ice_floes[1].fixed = true
-@test 1 == length(sim.overlaps)
@test 1 == length(sim.contact_pairs)
@test_approx_eq [1, 2] sim.contact_pairs[1]
-@test_approx_eq [2., 0.] sim.overlaps[1]
info("Testing findContacts(...)")
sim = deepcopy(sim_copy)
SeaIce.findContacts!(sim)
-@test 1 == length(sim.overlaps)
@test 1 == length(sim.contact_pairs)
@test_approx_eq [1, 2] sim.contact_pairs[1]
-@test_approx_eq [2., 0.] sim.overlaps[1]
@test_throws ErrorException SeaIce.findContacts!(sim, method="")
@@ -51,20 +45,17 @@ sim = deepcopy(sim_copy)
sim.ice_floes[1].fixed = true
sim.ice_floes[2].fixed = true
SeaIce.findContacts!(sim)
-@test 0 == length(sim.overlaps)
@test 0 == length(sim.contact_pairs)
sim = deepcopy(sim_copy)
SeaIce.disableIceFloe!(sim, 1)
SeaIce.findContacts!(sim)
-@test 0 == length(sim.overlaps)
@test 0 == length(sim.contact_pairs)
sim = deepcopy(sim_copy)
SeaIce.disableIceFloe!(sim, 1)
SeaIce.disableIceFloe!(sim, 2)
SeaIce.findContacts!(sim)
-@test 0 == length(sim.overlaps)
@test 0 == length(sim.contact_pairs)
info("Testing if interact(...) removes contacts correctly")
@@ -73,10 +64,8 @@ SeaIce.findContacts!(sim)
SeaIce.interact!(sim)
SeaIce.findContacts!(sim)
-@test 1 == length(sim.overlaps)
@test 1 == length(sim.contact_pairs)
@test_approx_eq [1, 2] sim.contact_pairs[1]
-@test_approx_eq [2., 0.] sim.overlaps[1]
info("Testing findContactsOceanGrid(...)")
sim = deepcopy(sim_copy)
@@ -84,10 +73,8 @@ sim.ocean = SeaIce.createRegularOceanGrid([4, 4, 2], [80., 80., 2.])
SeaIce.sortIceFloesInOceanGrid!(sim)
SeaIce.findContactsOceanGrid!(sim)
-@test 1 == length(sim.overlaps)
@test 1 == length(sim.contact_pairs)
@test_approx_eq [1, 2] sim.contact_pairs[1]
-@test_approx_eq [2., 0.] sim.overlaps[1]
sim = deepcopy(sim_copy)
sim.ocean = SeaIce.createRegularOceanGrid([4, 4, 2], [80., 80., 2.])
@@ -95,10 +82,8 @@ sim.ice_floes[1].fixed = true
SeaIce.sortIceFloesInOceanGrid!(sim)
SeaIce.findContactsOceanGrid!(sim)
-@test 1 == length(sim.overlaps)
@test 1 == length(sim.contact_pairs)
@test_approx_eq [1, 2] sim.contact_pairs[1]
-@test_approx_eq [2., 0.] sim.overlaps[1]
sim = deepcopy(sim_copy)
sim.ocean = SeaIce.createRegularOceanGrid([4, 4, 2], [80., 80., 2.])
@@ -107,7 +92,6 @@ sim.ice_floes[2].fixed = true
SeaIce.sortIceFloesInOceanGrid!(sim)
SeaIce.findContactsOceanGrid!(sim)
-@test 0 == length(sim.overlaps)
@test 0 == length(sim.contact_pairs)
info("Testing findContacts(...)")
@@ -116,14 +100,10 @@ sim.ocean = SeaIce.createRegularOceanGrid([4, 4, 2], [80., 80., 2.])
SeaIce.sortIceFloesInOceanGrid!(sim)
SeaIce.findContacts!(sim)
-@test 1 == length(sim.overlaps)
@test 1 == length(sim.contact_pairs)
@test_approx_eq [1, 2] sim.contact_pairs[1]
-@test_approx_eq [2., 0.] sim.overlaps[1]
-@test 1 == length(sim.overlaps)
@test 1 == length(sim.contact_pairs)
@test_approx_eq [1, 2] sim.contact_pairs[1]
-@test_approx_eq [2., 0.] sim.overlaps[1]
@test_throws ErrorException SeaIce.findContacts!(sim, method="")
