tUse Young's modulus and Poisson's ratio by default for contact mechanics - 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 7fe5b27334a798544f5a5ae76b8c259026d37e65
(DIR) parent b4856096f45707589baa494ad735de2128f7ec35
(HTM) Author: Anders Damsgaard <andersd@riseup.net>
Date: Wed, 24 May 2017 15:09:43 -0400
Use Young's modulus and Poisson's ratio by default for contact mechanics
Diffstat:
M src/icefloe.jl | 5 +++--
M src/interaction.jl | 70 ++++++++++++++++++++-----------
M src/temporal.jl | 19 +++++++++++++++----
M test/collision-2floes-oblique.jl | 22 ++++++++++++----------
M test/vtk.jl | 2 +-
5 files changed, 77 insertions(+), 41 deletions(-)
---
(DIR) diff --git a/src/icefloe.jl b/src/icefloe.jl
t@@ -21,13 +21,14 @@ function addIceFloeCylindrical(simulation::Simulation,
ang_acc::float = 0.,
torque::float = 0.,
density::float = 934.,
- contact_stiffness_normal::float = 1.e6,
+ contact_stiffness_normal::float = 1e7,
contact_stiffness_tangential::float = 0.,
contact_viscosity_normal::float = 0.,
contact_viscosity_tangential::float = 0.,
contact_static_friction::float = 0.4,
contact_dynamic_friction::float = 0.4,
- youngs_modulus::float = 2e9, # Hopkins 2004
+ youngs_modulus::float = 2e7,
+ #youngs_modulus::float = 2e9, # Hopkins 2004
poissons_ratio::float = 0.185, # Hopkins 2004
tensile_strength::float = 500e3, # Hopkins 2004
compressive_strength_prefactor::float = 1285e3,
(DIR) diff --git a/src/interaction.jl b/src/interaction.jl
t@@ -33,9 +33,14 @@ end
export interactIceFloes!
"""
+ interactIceFloes!(simulation::Simulation, i::Int, j::Int, ic::Int)
+
Resolve an grain-to-grain interaction using a prescibed contact law. This
function adds the compressive force of the interaction to the ice floe
`pressure` field of mean compressive stress on the ice floe sides.
+
+The function uses the macroscopic contact parameterization based on Young's
+modulus and Poisson's ratio if Young's modulus is a positive value.
"""
function interactIceFloes!(simulation::Simulation, i::Int, j::Int, ic::Int)
t@@ -80,19 +85,37 @@ function interactIceFloes!(simulation::Simulation, i::Int, j::Int, ic::Int)
abs(delta_n)/2.
# Contact mechanical parameters
- k_n_harmonic_mean =
- harmonicMean(simulation.ice_floes[i].contact_stiffness_normal,
- simulation.ice_floes[j].contact_stiffness_normal)
-
- k_t_harmonic_mean =
- harmonicMean(simulation.ice_floes[i].contact_stiffness_tangential,
- simulation.ice_floes[j].contact_stiffness_tangential)
+ if simulation.ice_floes[i].youngs_modulus > 0. &&
+ simulation.ice_floes[j].youngs_modulus > 0.
+
+ E = harmonicMean(simulation.ice_floes[i].youngs_modulus,
+ simulation.ice_floes[j].youngs_modulus)
+ ν = 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
+ k_t = k_n*2.*(1. - ν^2.)/((2. - ν)*(1. + ν)) # Obermayr et al 2011
+
+ else # Micromechanical parameterization: k_n and k_t set explicitly
+ k_n = harmonicMean(simulation.ice_floes[i].contact_stiffness_normal,
+ simulation.ice_floes[j].contact_stiffness_normal)
+
+ k_t =
+ harmonicMean(simulation.ice_floes[i].contact_stiffness_tangential,
+ simulation.ice_floes[j].contact_stiffness_tangential)
+ end
- gamma_n_harmonic_mean = harmonicMean(
+ gamma_n = harmonicMean(
simulation.ice_floes[i].contact_viscosity_normal,
simulation.ice_floes[j].contact_viscosity_normal)
- gamma_t_harmonic_mean = harmonicMean(
+ gamma_t = harmonicMean(
simulation.ice_floes[i].contact_viscosity_tangential,
simulation.ice_floes[j].contact_viscosity_tangential)
t@@ -100,28 +123,28 @@ function interactIceFloes!(simulation::Simulation, i::Int, j::Int, ic::Int)
simulation.ice_floes[j].contact_dynamic_friction)
# Determine contact forces
- if k_n_harmonic_mean ≈ 0. && gamma_n_harmonic_mean ≈ 0.
+ if k_n ≈ 0. && gamma_n ≈ 0.
force_n = 0.
- elseif k_n_harmonic_mean > 0. && gamma_n_harmonic_mean ≈ 0.
- force_n = -k_n_harmonic_mean*delta_n
+ elseif k_n > 0. && gamma_n ≈ 0.
+ force_n = -k_n*delta_n
- elseif k_n_harmonic_mean > 0. && gamma_n_harmonic_mean > 0.
- force_n = -k_n_harmonic_mean*delta_n - gamma_n_harmonic_mean*vel_n
+ elseif k_n > 0. && gamma_n > 0.
+ force_n = -k_n*delta_n - gamma_n*vel_n
if force_n < 0.
force_n = 0.
end
else
- error("unknown contact_normal_rheology (k_n = $k_n_harmonic_mean," *
- " gamma_n = $gamma_n_harmonic_mean")
+ error("unknown contact_normal_rheology (k_n = $k_n," *
+ " gamma_n = $gamma_n")
end
- if k_t_harmonic_mean ≈ 0. && gamma_t_harmonic_mean ≈ 0.
+ if k_t ≈ 0. && gamma_t ≈ 0.
# do nothing
- elseif k_t_harmonic_mean ≈ 0. && gamma_t_harmonic_mean > 0.
- force_t = abs(gamma_t_harmonic_mean * vel_t)
+ elseif k_t ≈ 0. && gamma_t > 0.
+ force_t = abs(gamma_t * vel_t)
if force_t > mu_d_minimum*abs(force_n)
force_t = mu_d_minimum*abs(force_n)
end
t@@ -129,19 +152,18 @@ function interactIceFloes!(simulation::Simulation, i::Int, j::Int, ic::Int)
force_t = -force_t
end
- elseif k_t_harmonic_mean > 0.
+ elseif k_t > 0.
- force_t = -k_t_harmonic_mean*delta_t - gamma_t_harmonic_mean*vel_t
+ force_t = -k_t*delta_t - gamma_t*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
+ delta_t = (-force_t - gamma_t*vel_t)/k_t
end
else
error("unknown contact_tangential_rheology (k_t = " *
- "$k_t_harmonic_mean, gamma_t = $gamma_t_harmonic_mean")
+ "$k_t, gamma_t = $gamma_t")
end
end
simulation.ice_floes[i].contact_parallel_displacement[ic] = delta_t*t
(DIR) diff --git a/src/temporal.jl b/src/temporal.jl
t@@ -99,13 +99,24 @@ function findLargestIceFloeStiffness(simulation::Simulation)
i_n_max = -1
i_t_max = -1
for i in length(simulation.ice_floes)
+
icefloe = simulation.ice_floes[i]
- if icefloe.contact_stiffness_normal > k_n_max
- k_n_max = icefloe.contact_stiffness_normal
+
+ if icefloe.youngs_modulus > 0.
+ k_n = icefloe.youngs_modulus*icefloe.thickness # A = h*r
+ k_t = k_n*2.*(1. - icefloe.poissons_ratio^2.)/
+ ((2. - icefloe.poissons_ratio)*(1. + icefloe.poissons_ratio))
+ else
+ k_n = icefloe.contact_stiffness_normal
+ k_t = icefloe.contact_stiffness_tangential
+ end
+
+ if k_n > k_n_max
+ k_n_max = k_n
i_n_max = i
end
- if icefloe.contact_stiffness_tangential > k_t_max
- k_t_max = icefloe.contact_stiffness_tangential
+ if k_t > k_t_max
+ k_t_max = k_t
i_t_max = i
end
end
(DIR) diff --git a/test/collision-2floes-oblique.jl b/test/collision-2floes-oblique.jl
t@@ -13,6 +13,8 @@ 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[1].contact_dynamic_friction = 0.
+sim.ice_floes[2].contact_dynamic_friction = 0.
sim.ice_floes[2].fixed = true
E_kin_lin_init = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
t@@ -68,6 +70,8 @@ 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
+sim.ice_floes[1].contact_dynamic_friction = 0.
+sim.ice_floes[2].contact_dynamic_friction = 0.
E_kin_lin_init = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
E_kin_rot_init = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
t@@ -229,7 +233,7 @@ E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
info("Testing kinetic energy conservation with Two-term Taylor scheme")
sim = deepcopy(sim_init)
SeaIce.setTimeStep!(sim, epsilon=0.007)
-tol = 0.02
+tol = 0.04
info("Relative tolerance: $(tol*100.)%")
SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
verbose=verbose)
t@@ -242,7 +246,7 @@ E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
info("Testing kinetic energy conservation with Three-term Taylor scheme")
sim = deepcopy(sim_init)
SeaIce.setTimeStep!(sim, epsilon=0.07)
-tol = 0.02
+tol = 0.04
info("Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)")
SeaIce.run!(sim, temporal_integration_method="Three-term Taylor",
verbose=verbose)
t@@ -292,7 +296,7 @@ E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
info("Testing kinetic energy conservation with Two-term Taylor scheme")
sim = deepcopy(sim_init)
SeaIce.setTimeStep!(sim, epsilon=0.007)
-tol = 0.02
+tol = 0.04
info("Relative tolerance: $(tol*100.)%")
SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
verbose=verbose)
t@@ -305,7 +309,7 @@ E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
info("Testing kinetic energy conservation with Three-term Taylor scheme")
sim = deepcopy(sim_init)
SeaIce.setTimeStep!(sim, epsilon=0.07)
-tol = 0.02
+tol = 0.04
info("Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)")
SeaIce.run!(sim, temporal_integration_method="Three-term Taylor",
verbose=verbose)
t@@ -325,8 +329,6 @@ sim = SeaIce.createSimulation(id="test")
SeaIce.addIceFloeCylindrical(sim, [0., 0.], 10., 1., verbose=verbose)
SeaIce.addIceFloeCylindrical(sim, [19.0, -10.], 10., 1., verbose=verbose)
sim.ice_floes[2].lin_vel[1] = -0.1
-sim.ice_floes[1].contact_viscosity_tangential = 1e4
-sim.ice_floes[2].contact_viscosity_tangential = 1e4
E_kin_lin_init = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
E_kin_rot_init = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
t@@ -355,7 +357,7 @@ E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
info("Testing kinetic energy conservation with Two-term Taylor scheme")
sim = deepcopy(sim_init)
SeaIce.setTimeStep!(sim, epsilon=0.007)
-tol = 0.02
+tol = 0.04
info("Relative tolerance: $(tol*100.)%")
SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
verbose=verbose)
t@@ -368,7 +370,7 @@ E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
info("Testing kinetic energy conservation with Three-term Taylor scheme")
sim = deepcopy(sim_init)
SeaIce.setTimeStep!(sim, epsilon=0.07)
-tol = 0.02
+tol = 0.04
info("Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)")
SeaIce.run!(sim, temporal_integration_method="Three-term Taylor",
verbose=verbose)
t@@ -424,7 +426,7 @@ E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
info("Testing kinetic energy conservation with Two-term Taylor scheme")
sim = deepcopy(sim_init)
SeaIce.setTimeStep!(sim, epsilon=0.007)
-tol = 0.02
+tol = 0.04
info("Relative tolerance: $(tol*100.)%")
SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
verbose=verbose)
t@@ -437,7 +439,7 @@ E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
info("Testing kinetic energy conservation with Three-term Taylor scheme")
sim = deepcopy(sim_init)
SeaIce.setTimeStep!(sim, epsilon=0.07)
-tol = 0.03
+tol = 0.04
info("Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)")
SeaIce.run!(sim, temporal_integration_method="Three-term Taylor",
verbose=verbose)
(DIR) diff --git a/test/vtk.jl b/test/vtk.jl
t@@ -26,7 +26,7 @@ else
end
icefloechecksum =
-"4885bc7c0eccc5e43a77d370295841fae710465c8bcc120a13681a0947ffbd53 " *
+"de1ee347a6a15c665d3fa122c73092b3aa200bf8180007303f43f71396094d8d " *
"test.icefloes.1.vtu\n"
oceanchecksum =