commit 8c6f1f493f21bb6512ce390951af2aa1e3d33bfd
parent 6e04488619808a281e658e6ca8a104354eacc98d
Author: Anders Damsgaard <andersd@riseup.net>
Date: Tue, 2 May 2017 14:03:11 -0400
implement tangential viscosity and friction
Diffstat:
4 files changed, 338 insertions(+), 37 deletions(-)
diff --git a/src/interaction.jl b/src/interaction.jl
@@ -5,7 +5,8 @@ export interact!
Resolve mechanical interaction between all particle pairs.
"""
function interact!(simulation::Simulation;
- contact_model::String = "Linear Elastic No Tangential")
+ contact_normal_rheology::String = "Linear Elastic",
+ contact_tangential_rheology::String = "None")
# IceFloe to grain collisions
while !isempty(simulation.contact_pairs)
@@ -15,7 +16,9 @@ function interact!(simulation::Simulation;
pop!(simulation.contact_parallel_displacement)
interactIceFloes!(simulation, contact_pair[1], contact_pair[2],
overlap_vector, contact_parallel_displacement,
- contact_model=contact_model)
+ contact_normal_rheology=contact_normal_rheology,
+ contact_tangential_rheology=
+ contact_tangential_rheology)
end
end
@@ -26,34 +29,42 @@ function adds the compressive force of the interaction to the ice floe
`pressure` field of mean compressive stress on the ice floe sides.
"""
function interactIceFloes!(simulation::Simulation,
- i::Integer, j::Integer,
+ i::Int, j::Int,
overlap_vector::Array{Float64, 1},
contact_parallel_displacement::Array{Float64, 1};
- contact_model::String =
- "Linear Elastic No Tangential")
+ contact_normal_rheology::String = "Linear Elastic",
+ contact_tangential_rheology::String = "None")
- force = zeros(2)
+ force_n = zeros(2)
+ force_t = zeros(2)
- if contact_model == "None"
- return nothing
-
- elseif contact_model == "Linear Elastic No Tangential" ||
- contact_model == "LinearElastic"
+ if contact_normal_rheology == "None"
+ # do nothing
+ elseif contact_normal_rheology == "Linear Elastic"
force_n = interactNormalLinearElastic(simulation, i, j, overlap_vector)
+ else
+ error("Unknown contact_normal_rheology '$contact_normal_rheology'")
+ end
- if contact_model == "Linear Elastic"
- force_t, torque = interactTangentialLinearElasticFrictional(
- simulation, i, j,
- overlap_vector,
- contact_parallel_displacement)
- end
-
+ if contact_tangential_rheology == "None"
+ # do nothing
+ elseif contact_tangential_rheology == "Linear Viscous Frictional"
+ force_t = interactTangentialLinearViscousFrictional(simulation, i, j,
+ overlap_vector,
+ force_n)
else
- error("Unknown contact_model interaction model '$contact_model'")
+ error("Unknown contact_tangential_rheology ",
+ "'$contact_tangential_rheology'")
end
- simulation.ice_floes[i].force += force_n;
- simulation.ice_floes[j].force -= force_n;
+ simulation.ice_floes[i].force += force_n + force_t;
+ simulation.ice_floes[j].force -= force_n + force_t;
+
+ if norm(force_t) > 0.
+ torque = -findTorque(simulation, overlap_vector, force_t, i, j)
+ simulation.ice_floes[i].torque += torque
+ simulation.ice_floes[j].torque += torque
+ end
simulation.ice_floes[i].pressure +=
norm(force_n)/simulation.ice_floes[i].side_surface_area;
@@ -67,7 +78,7 @@ Resolves linear-elastic interaction between two ice floes in the contact-normal
direction.
"""
function interactNormalLinearElastic(simulation::Simulation,
- i::Integer, j::Integer,
+ i::Int, j::Int,
overlap_vector::vector)
k_n_harmonic_mean =
@@ -82,25 +93,56 @@ export interactTangentialLinearElastic
Resolves linear-elastic interaction between two ice floes in the
contact-parallel (tangential) direction.
"""
-function interactTangentialLinearElasticFrictional(simulation::Simulation,
- i::Integer, j::Integer,
- overlap_vector::vector,
- contact_parallel_displacement::vector)
+function interactTangentialLinearViscousFrictional(simulation::Simulation,
+ i::Int, j::Integer,
+ overlap_vector::vector,
+ force_n::vector)
- k_t_harmonic_mean =
- harmonicMean(simulation.ice_floes[i].contact_stiffness_tangential,
- simulation.ice_floes[j].contact_stiffness_tangential)
+ contact_parallel_velocity = findContactParallelVelocity(simulation, i, j,
+ overlap_vector)
- # correct displacement for contact rotation
- n = overlap_vector/norm(overlap_vector)
- contact_parallel_displacement -= (n * dot(n, contact_parallel_displacement))
+ if norm(contact_parallel_velocity) ≈ 0.
+ return [0., 0.]
+ end
+
+ gamma_t_harmonic_mean = harmonicMean(
+ simulation.ice_floes[i].contact_viscosity_tangential,
+ simulation.ice_floes[j].contact_viscosity_tangential)
- force_t = k_t_harmonic_mean * contact_parallel_displacement
+ if norm(gamma_t_harmonic_mean) ≈ 0.
+ return [0., 0.]
+ end
- return force, torque, contact_parallel_displacement
+ mu_d_minimum = min(simulation.ice_floes[i].contact_dynamic_friction,
+ simulation.ice_floes[j].contact_dynamic_friction)
+ return -min(norm(gamma_t_harmonic_mean*contact_parallel_velocity),
+ mu_d_minimum*norm(force_n))*
+ contact_parallel_velocity/norm(contact_parallel_velocity)
end
function harmonicMean(a::Any, b::Any)
return 2.*a*b/(a + b)
end
+
+function findTorque(simulation::Simulation, overlap_vector::vector,
+ force_t::vector, i::Int, j::Int)
+ n = overlap_vector/norm(overlap_vector)
+ return -findEffectiveRadius(simulation, i, j, overlap_vector)*
+ (n[1]*force_t[2] - n[2]*force_t[1])
+end
+
+function findEffectiveRadius(simulation::Simulation, i::Int, j::Int,
+ overlap_vector::vector)
+ return harmonicMean(simulation.ice_floes[i].contact_radius,
+ simulation.ice_floes[j].contact_radius)
+ - norm(overlap_vector)/2.
+end
+
+function findContactParallelVelocity(simulation::Simulation, i::Int, j::Int,
+ overlap_vector::vector)
+ return simulation.ice_floes[i].lin_vel -
+ simulation.ice_floes[j].lin_vel +
+ findEffectiveRadius(simulation, i, j, overlap_vector)*
+ (simulation.ice_floes[i].ang_vel + simulation.ice_floes[j].ang_vel)
+end
diff --git a/src/simulation.jl b/src/simulation.jl
@@ -89,8 +89,8 @@ function run!(simulation::Simulation;
show_file_output::Bool=true,
single_step::Bool=false,
temporal_integration_method::String="Three-term Taylor",
- contact_model::String="Three-term Taylor",
- )
+ contact_normal_rheology::String = "Linear Elastic",
+ contact_tangential_rheology::String = "None")
if single_step && simulation.time >= simulation.time_total
simulation.time_total += simulation.time_step
@@ -122,7 +122,9 @@ function run!(simulation::Simulation;
else
findContacts!(simulation, method="all to all")
end
- interact!(simulation)
+ interact!(simulation,
+ contact_normal_rheology=contact_normal_rheology,
+ contact_tangential_rheology=contact_tangential_rheology)
if typeof(simulation.ocean.input_file) != Bool
addOceanDrag!(simulation)
end
diff --git a/test/collision-2floes-oblique.jl b/test/collision-2floes-oblique.jl
@@ -0,0 +1,256 @@
+#!/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
+
+info("## Contact-normal elasticity only")
+info("# One ice floe fixed")
+sim = SeaIce.createSimulation(id="test")
+SeaIce.addIceFloeCylindrical(sim, [0., 10.], 10., 1., verbose=verbose)
+SeaIce.addIceFloeCylindrical(sim, [19., 0.], 10., 1., verbose=verbose)
+sim.ice_floes[1].lin_vel[1] = 0.1
+sim.ice_floes[2].fixed = true
+
+E_kin_lin_init = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_init = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+
+# With decreasing timestep (epsilon towards 0), the explicit integration scheme
+# should become more correct
+
+SeaIce.setTotalTime!(sim, 30.0)
+#sim.file_time_step = 1.
+sim_init = deepcopy(sim)
+
+info("Testing kinetic energy conservation with Two-term Taylor scheme")
+SeaIce.setTimeStep!(sim, epsilon=0.07)
+tol = 0.2
+info("Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)")
+SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
+ contact_tangential_rheology="None", verbose=verbose)
+
+E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+@test_approx_eq_eps E_kin_lin_init E_kin_lin_final E_kin_lin_init*tol
+@test_approx_eq E_kin_rot_init E_kin_rot_final
+
+
+info("Testing kinetic energy conservation with Two-term Taylor scheme")
+sim = deepcopy(sim_init)
+SeaIce.setTimeStep!(sim, epsilon=0.007)
+tol = 0.02
+info("Relative tolerance: $(tol*100.)%")
+SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
+ contact_tangential_rheology="None", verbose=verbose)
+
+E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+@test_approx_eq_eps E_kin_lin_init E_kin_lin_final E_kin_lin_init*tol
+@test_approx_eq E_kin_rot_init E_kin_rot_final
+
+
+info("Testing kinetic energy conservation with Three-term Taylor scheme")
+sim = deepcopy(sim_init)
+SeaIce.setTimeStep!(sim, epsilon=0.07)
+tol = 0.01
+info("Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)")
+SeaIce.run!(sim, temporal_integration_method="Three-term Taylor",
+ contact_tangential_rheology="None", verbose=verbose)
+
+E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+@test_approx_eq_eps E_kin_lin_init E_kin_lin_final E_kin_lin_init*tol
+@test_approx_eq E_kin_rot_init E_kin_rot_final
+
+info("# Ice floes free to move")
+
+sim = SeaIce.createSimulation(id="test")
+SeaIce.addIceFloeCylindrical(sim, [0., 10.], 10., 1., verbose=verbose)
+SeaIce.addIceFloeCylindrical(sim, [19.0, 0.], 10., 1., verbose=verbose)
+sim.ice_floes[1].lin_vel[1] = 0.1
+
+E_kin_lin_init = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_init = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+
+# With decreasing timestep (epsilon towards 0), the explicit integration scheme
+# should become more correct
+
+SeaIce.setTotalTime!(sim, 30.0)
+sim_init = deepcopy(sim)
+
+info("Testing kinetic energy conservation with Two-term Taylor scheme")
+SeaIce.setTimeStep!(sim, epsilon=0.07)
+tol = 0.2
+info("Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)")
+SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
+ contact_tangential_rheology="None", verbose=verbose)
+
+E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+@test_approx_eq_eps E_kin_lin_init E_kin_lin_final E_kin_lin_init*tol
+@test_approx_eq E_kin_rot_init E_kin_rot_final
+
+
+info("Testing kinetic energy conservation with Two-term Taylor scheme")
+sim = deepcopy(sim_init)
+SeaIce.setTimeStep!(sim, epsilon=0.007)
+tol = 0.02
+info("Relative tolerance: $(tol*100.)%")
+SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
+ contact_tangential_rheology="None", verbose=verbose)
+
+E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+@test_approx_eq_eps E_kin_lin_init E_kin_lin_final E_kin_lin_init*tol
+@test_approx_eq E_kin_rot_init E_kin_rot_final
+
+
+info("Testing kinetic energy conservation with Three-term Taylor scheme")
+sim = deepcopy(sim_init)
+SeaIce.setTimeStep!(sim, epsilon=0.07)
+tol = 0.01
+info("Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)")
+SeaIce.run!(sim, temporal_integration_method="Three-term Taylor",
+ contact_tangential_rheology="None", verbose=verbose)
+
+E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+@test_approx_eq_eps E_kin_lin_init E_kin_lin_final E_kin_lin_init*tol
+@test_approx_eq E_kin_rot_init E_kin_rot_final
+
+
+info("## Contact-normal elasticity and tangential viscosity and friction")
+info("# One ice floe fixed")
+sim = SeaIce.createSimulation(id="test")
+SeaIce.addIceFloeCylindrical(sim, [0., 10.], 10., 1., verbose=verbose)
+SeaIce.addIceFloeCylindrical(sim, [19., 0.], 10., 1., verbose=verbose)
+sim.ice_floes[1].lin_vel[1] = 0.1
+sim.ice_floes[2].fixed = true
+
+E_kin_lin_init = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_init = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+
+# With decreasing timestep (epsilon towards 0), the explicit integration scheme
+# should become more correct
+
+SeaIce.setTotalTime!(sim, 30.0)
+#sim.file_time_step = 1.
+sim_init = deepcopy(sim)
+
+info("Testing kinetic energy conservation with Two-term Taylor scheme")
+SeaIce.setTimeStep!(sim, epsilon=0.07)
+tol = 0.5
+info("Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)")
+SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
+ contact_tangential_rheology="Linear Viscous Frictional",
+ verbose=verbose)
+
+@test sim.ice_floes[1].ang_pos > 0.
+@test sim.ice_floes[1].ang_vel > 0.
+@test sim.ice_floes[2].ang_pos ≈ 0.
+@test sim.ice_floes[2].ang_vel ≈ 0.
+E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+@test_approx_eq_eps E_kin_lin_init+E_kin_rot_init E_kin_lin_final+E_kin_rot_final E_kin_lin_init*tol
+
+
+info("Testing kinetic energy conservation with Two-term Taylor scheme")
+sim = deepcopy(sim_init)
+SeaIce.setTimeStep!(sim, epsilon=0.007)
+tol = 0.2
+info("Relative tolerance: $(tol*100.)%")
+SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
+ contact_tangential_rheology="Linear Viscous Frictional",
+ verbose=verbose)
+
+@test sim.ice_floes[1].ang_pos > 0.
+@test sim.ice_floes[1].ang_vel > 0.
+@test sim.ice_floes[2].ang_pos ≈ 0.
+@test sim.ice_floes[2].ang_vel ≈ 0.
+E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+@test_approx_eq_eps E_kin_lin_init+E_kin_rot_init E_kin_lin_final+E_kin_rot_final E_kin_lin_init*tol
+
+
+info("Testing kinetic energy conservation with Three-term Taylor scheme")
+sim = deepcopy(sim_init)
+SeaIce.setTimeStep!(sim, epsilon=0.07)
+tol = 0.2
+info("Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)")
+SeaIce.run!(sim, temporal_integration_method="Three-term Taylor",
+ contact_tangential_rheology="Linear Viscous Frictional",
+ verbose=verbose)
+
+@test sim.ice_floes[1].ang_pos > 0.
+@test sim.ice_floes[1].ang_vel > 0.
+@test sim.ice_floes[2].ang_pos ≈ 0.
+@test sim.ice_floes[2].ang_vel ≈ 0.
+E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+@test_approx_eq_eps E_kin_lin_init+E_kin_rot_init E_kin_lin_final+E_kin_rot_final E_kin_lin_init*tol
+
+info("# Ice floes free to move")
+
+sim = SeaIce.createSimulation(id="test")
+SeaIce.addIceFloeCylindrical(sim, [0., 10.], 10., 1., verbose=verbose)
+SeaIce.addIceFloeCylindrical(sim, [19.0, 0.], 10., 1., verbose=verbose)
+sim.ice_floes[1].lin_vel[1] = 0.1
+
+E_kin_lin_init = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_init = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+
+# With decreasing timestep (epsilon towards 0), the explicit integration scheme
+# should become more correct
+
+SeaIce.setTotalTime!(sim, 30.0)
+sim_init = deepcopy(sim)
+
+info("Testing kinetic energy conservation with Two-term Taylor scheme")
+SeaIce.setTimeStep!(sim, epsilon=0.07)
+tol = 0.2
+info("Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)")
+SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
+ contact_tangential_rheology="Linear Viscous Frictional",
+ verbose=verbose)
+
+@test sim.ice_floes[1].ang_pos > 0.
+@test sim.ice_floes[1].ang_vel > 0.
+@test sim.ice_floes[2].ang_pos > 0.
+@test sim.ice_floes[2].ang_vel > 0.
+E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+@test_approx_eq_eps E_kin_lin_init+E_kin_rot_init E_kin_lin_final+E_kin_rot_final E_kin_lin_init*tol
+
+info("Testing kinetic energy conservation with Two-term Taylor scheme")
+sim = deepcopy(sim_init)
+SeaIce.setTimeStep!(sim, epsilon=0.007)
+tol = 0.1
+info("Relative tolerance: $(tol*100.)%")
+SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
+ contact_tangential_rheology="Linear Viscous Frictional",
+ verbose=verbose)
+
+E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+@test_approx_eq_eps E_kin_lin_init+E_kin_rot_init E_kin_lin_final+E_kin_rot_final E_kin_lin_init*tol
+
+
+info("Testing kinetic energy conservation with Three-term Taylor scheme")
+sim = deepcopy(sim_init)
+SeaIce.setTimeStep!(sim, epsilon=0.07)
+tol = 0.1
+info("Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)")
+SeaIce.run!(sim, temporal_integration_method="Three-term Taylor",
+ contact_tangential_rheology="Linear Viscous Frictional",
+ verbose=verbose)
+
+@test sim.ice_floes[1].ang_pos > 0.
+@test sim.ice_floes[1].ang_vel > 0.
+@test sim.ice_floes[2].ang_pos > 0.
+@test sim.ice_floes[2].ang_vel > 0.
+E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+@test_approx_eq_eps E_kin_lin_init+E_kin_rot_init E_kin_lin_final+E_kin_rot_final E_kin_lin_init*tol
diff --git a/test/runtests.jl b/test/runtests.jl
@@ -3,6 +3,7 @@ using Base.Test
include("contact-search-and-geometry.jl")
include("collision-2floes-normal.jl")
+include("collision-2floes-oblique.jl")
include("netcdf.jl")
include("vtk.jl")
include("grid.jl")
