commit 6250a89c16baa8c912f836f2586cc9dc8e9537b6
parent 2e88c8f814b4838ec6ec179feb6a565d46a7b155
Author: Anders Damsgaard <andersd@riseup.net>
Date: Thu, 25 May 2017 14:03:42 -0400
linearly increase cohesion with contact time until hard limit, remove contacts within interaction function
Diffstat:
6 files changed, 60 insertions(+), 12 deletions(-)
diff --git a/src/contact_search.jl b/src/contact_search.jl
@@ -149,6 +149,7 @@ function checkAndAddContact!(sim::Simulation, i::Int, j::Int)
sim.ice_floes[i].contacts[ic] = j
sim.ice_floes[i].contact_parallel_displacement[ic] =
zeros(2)
+ sim.ice_floes[i].contact_age[ic] = 0.
break
end
end
diff --git a/src/datatypes.jl b/src/datatypes.jl
@@ -61,6 +61,7 @@ type IceFloeCylindrical
ocean_grid_pos::Array{Int, 1}
contacts::Array{Int, 1}
contact_parallel_displacement::Array{Array{Float64, 1}, 1}
+ contact_age::Array{Float64, 1}
end
# Type for gathering data from ice floe objects into single arrays
diff --git a/src/icefloe.jl b/src/icefloe.jl
@@ -46,9 +46,10 @@ function addIceFloeCylindrical(simulation::Simulation,
n_contacts::Int = 0,
contacts::Array{Int, 1} = zeros(Int, Nc_max),
contact_parallel_displacement::
- Array{Array{Float64, 1}, 1}
- =
- Array{Array{Float64, 1}, 1}(Nc_max))
+ Array{Array{Float64, 1}, 1} =
+ Array{Array{Float64, 1}, 1}(Nc_max),
+ contact_age::Array{Float64, 1} =
+ zeros(Float64, Nc_max))
# Check input values
if length(lin_pos) != 2
@@ -126,7 +127,8 @@ function addIceFloeCylindrical(simulation::Simulation,
n_contacts,
ocean_grid_pos,
contacts,
- contact_parallel_displacement
+ contact_parallel_displacement,
+ contact_age
)
# Overwrite previous placeholder values
diff --git a/src/interaction.jl b/src/interaction.jl
@@ -16,6 +16,7 @@ function interact!(simulation::Simulation)
continue
end
+ """
if norm(simulation.ice_floes[i].lin_pos -
simulation.ice_floes[j].lin_pos) -
(simulation.ice_floes[i].contact_radius +
@@ -25,8 +26,9 @@ function interact!(simulation::Simulation)
simulation.ice_floes[i].n_contacts -= 1
simulation.ice_floes[j].n_contacts -= 1
else
- interactIceFloes!(simulation, i, j, ic)
- end
+ """
+ interactIceFloes!(simulation, i, j, ic)
+ #end
end
end
end
@@ -76,6 +78,9 @@ function interactIceFloes!(simulation::Simulation, i::Int, j::Int, ic::Int)
# Effective radius
R_ij = harmonicMean(simulation.ice_floes[i].contact_radius,
simulation.ice_floes[j].contact_radius) - abs(δ_n)/2.
+ # Contact area
+ A_ij = R_ij*min(simulation.ice_floes[i].thickness,
+ simulation.ice_floes[j].thickness)
# Contact mechanical parameters
if simulation.ice_floes[i].youngs_modulus > 0. &&
@@ -86,10 +91,6 @@ function interactIceFloes!(simulation::Simulation, i::Int, j::Int, ic::Int)
ν = harmonicMean(simulation.ice_floes[i].poissons_ratio,
simulation.ice_floes[j].poissons_ratio)
- # Contact area
- A_ij = R_ij*min(simulation.ice_floes[i].thickness,
- simulation.ice_floes[j].thickness)
-
# Effective normal and tangential stiffness
k_n = E*A_ij/R_ij
#k_t = k_n*ν # Kneib et al 2016
@@ -126,8 +127,27 @@ function interactIceFloes!(simulation::Simulation, i::Int, j::Int, ic::Int)
end
else
- error("unknown contact_normal_rheology (k_n = $k_n," *
- " γ_n = $γ_n")
+ error("unknown contact_normal_rheology (k_n = $k_n, γ_n = $γ_n")
+ end
+
+ # Contact tensile strength increases linearly with contact age until tensile
+ # stress exceeds tensile strength
+ if δ_n > 0.
+
+ # linearly increase tensile strength with time until max. value
+ tensile_strength = min(simulation.ice_floes[i].contact_age[ic]/
+ (60.*60.*24.), 1.)*
+ simulation.ice_floes[i].tensile_strength*
+ simulation.ice_floes[i].thickness
+
+ # break bond
+ if abs(force_n) > tensile_strength*A_ij
+ force_n = 0.
+ force_t = 0.
+ simulation.ice_floes[i].contacts[ic] = 0 # remove contact
+ simulation.ice_floes[i].n_contacts -= 1
+ simulation.ice_floes[j].n_contacts -= 1
+ end
end
if k_t ≈ 0. && γ_t ≈ 0.
@@ -156,6 +176,7 @@ function interactIceFloes!(simulation::Simulation, i::Int, j::Int, ic::Int)
end
simulation.ice_floes[i].contact_parallel_displacement[ic] = δ_t*t
+ simulation.ice_floes[i].contact_age[ic] += simulation.time_step
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/test/cohesion.jl b/test/cohesion.jl
@@ -0,0 +1,22 @@
+#!/usr/bin/env julia
+
+# Check for conservation of kinetic energy (=momentum) during a normal collision
+# between two ice cylindrical ice floes
+
+info("#### $(basename(@__FILE__)) ####")
+
+verbose=false
+
+sim_init = SeaIce.createSimulation(id="test")
+SeaIce.addIceFloeCylindrical(sim_init, [0., 0.], 10., 1., verbose=verbose)
+SeaIce.addIceFloeCylindrical(sim_init, [18., 0.], 10., 1., verbose=verbose)
+sim_init.ice_floes[1].youngs_modulus = 1e-5 # repulsion is negligible
+sim_init.ice_floes[2].youngs_modulus = 1e-5 # repulsion is negligible
+SeaIce.setTimeStep!(sim_init, verbose=verbose)
+
+info("# Check contact age scheme")
+sim = deepcopy(sim_init)
+SeaIce.setTotalTime!(sim, 10.)
+sim.time_step = 1.
+SeaIce.run!(sim, verbose=verbose)
+@test sim.ice_floes[1].contact_age[1] ≈ sim.time
diff --git a/test/runtests.jl b/test/runtests.jl
@@ -7,6 +7,7 @@ include("contact-search-and-geometry.jl")
include("collision-2floes-normal.jl")
include("collision-5floes-normal.jl")
include("collision-2floes-oblique.jl")
+include("cohesion.jl")
include("netcdf.jl")
include("vtk.jl")
include("grid.jl")
