commit b320d55a80891c9261ff58288e9c751826b581fa
parent 5ce691186cd1e035b45adafe40e4ba9067902886
Author: Anders Damsgaard <andersd@riseup.net>
Date: Tue, 19 Dec 2017 10:22:57 -0500
Add wall-normal viscosity in parallel to elasticity
Diffstat:
4 files changed, 145 insertions(+), 8 deletions(-)
diff --git a/src/datatypes.jl b/src/datatypes.jl
@@ -83,6 +83,7 @@ mutable struct WallLinearFrictionless
vel::Float64 # Velocity (constant when bc == "normal stress")
acc::Float64 # Acceleration (zero when bc == "velocity")
force::Float64 # Sum of normal forces on wall
+ contact_viscosity_normal::Float64 # Wall-normal contact viscosity
end
# Type for gathering data from grain objects into single arrays
diff --git a/src/interaction.jl b/src/interaction.jl
@@ -49,6 +49,9 @@ function interactWalls!(sim::Simulation)
orientation::Float64 = 0.0
δ_n::Float64 = 0.0
k_n::Float64 = 0.0
+ γ_n::Float64 = 0.0
+ vel_n::Float64 = 0.0
+ force_n::Float64 = 0.0
for iw=1:length(sim.walls)
for i=1:length(sim.grains)
@@ -57,10 +60,12 @@ function interactWalls!(sim::Simulation)
sim.grains[i].lin_pos) - sim.walls[iw].pos)
# get overlap distance by projecting grain position onto wall-normal
- # vector
+ # vector. Overlap when δ_n < 0.0
δ_n = abs(dot(sim.walls[iw].normal, sim.grains[i].lin_pos) -
sim.walls[iw].pos) - sim.grains[i].contact_radius
+ vel_n = dot(sim.walls[iw].normal, sim.grains[i].lin_vel)
+
if δ_n < 0.
if sim.grains[i].youngs_modulus > 0.
k_n = sim.grains[i].youngs_modulus * sim.grains[i].thickness
@@ -68,9 +73,31 @@ function interactWalls!(sim::Simulation)
k_n = sim.grains[i].contact_stiffness_normal
end
- sim.grains[i].force += -k_n * δ_n .* sim.walls[iw].normal .*
+ γ_n = sim.walls[iw].contact_viscosity_normal
+
+ if k_n ≈ 0. && γ_n ≈ 0. # No interaction
+ force_n = 0.
+
+ elseif k_n > 0. && γ_n ≈ 0. # Elastic (Hookean)
+ force_n = k_n * δ_n
+
+ elseif k_n > 0. && γ_n > 0. # Elastic-viscous (Kelvin-Voigt)
+ force_n = k_n * δ_n + γ_n * vel_n
+
+ # Make sure that the visous component doesn't dominate over
+ # elasticity
+ if force_n > 0.
+ force_n = 0.
+ end
+
+ else
+ error("unknown contact_normal_rheology " *
+ "(k_n = $k_n, γ_n = $γ_n")
+ end
+
+ sim.grains[i].force += -force_n .* sim.walls[iw].normal .*
orientation
- sim.walls[iw].force += k_n * δ_n * orientation
+ sim.walls[iw].force += force_n * orientation
end
end
end
@@ -162,13 +189,13 @@ function interactGrains!(simulation::Simulation, i::Int, j::Int, ic::Int)
simulation.grains[j].contact_dynamic_friction)
# Determine contact forces
- if k_n ≈ 0. && γ_n ≈ 0.
+ if k_n ≈ 0. && γ_n ≈ 0. # No interaction
force_n = 0.
- elseif k_n > 0. && γ_n ≈ 0.
+ elseif k_n > 0. && γ_n ≈ 0. # Elastic (Hookean)
force_n = -k_n*δ_n
- elseif k_n > 0. && γ_n > 0.
+ elseif k_n > 0. && γ_n > 0. # Elastic-viscous (Kelvin-Voigt)
force_n = -k_n*δ_n - γ_n*vel_n
if force_n < 0.
force_n = 0.
diff --git a/src/wall.jl b/src/wall.jl
@@ -35,6 +35,10 @@ The only required arguments are
parameter. Units: [m/s]
* `force::Float64=0.`: sum of normal forces on the wall from interaction with
grains [N].
+* `contact_viscosity_normal::Float64=0.`: viscosity to apply in parallel to
+ elasticity in interactions between wall and particles [N/(m/s)]. When this
+ term is larger than zero, the wall-grain interaction acts like a sink of
+ kinetic energy.
* `verbose::Bool=true`: show verbose information during function call.
# Examples
@@ -74,6 +78,7 @@ function addWallLinearFrictionless!(simulation::Simulation,
vel::Float64 = 0.,
acc::Float64 = 0.,
force::Float64 = 0.,
+ contact_viscosity_normal = 0.,
verbose::Bool=true)
# Check input values
@@ -143,7 +148,8 @@ function addWallLinearFrictionless!(simulation::Simulation,
normal_stress,
vel,
acc,
- force)
+ force,
+ contact_viscosity_normal)
# Add to simulation object
addWall!(simulation, wall, verbose)
diff --git a/test/wall.jl b/test/wall.jl
@@ -50,7 +50,7 @@ Granular.addWallLinearFrictionless!(sim, [0., 1.], 1., verbose=false)
@test Granular.getWallSurfaceArea(sim, 1) ≈ 20.0*2.0
@test Granular.getWallSurfaceArea(sim, 2) ≈ 10.0*2.0
-info("# Test wall-grain interaction")
+info("# Test wall-grain interaction: elastic")
info("Wall present but no contact")
sim = Granular.createSimulation()
@@ -118,6 +118,109 @@ Granular.interactWalls!(sim)
@test sim.grains[1].force[1] ≈ 0.
@test sim.grains[1].force[2] < 0.
+info("# Test wall-grain interaction: elastic-viscous")
+
+info("Wall present but no contact")
+sim = Granular.createSimulation()
+sim.ocean = Granular.createRegularOceanGrid([1, 1, 1], [10., 20., 1.0])
+Granular.addGrainCylindrical!(sim, [ 0., 0.], 1., 2., verbose=false)
+Granular.addWallLinearFrictionless!(sim, [1., 0.], -1.01, verbose=false)
+sim.walls[1].contact_viscosity_normal = 1e3
+Granular.setTimeStep!(sim, verbose=false)
+Granular.interactWalls!(sim)
+@test sim.walls[1].force ≈ 0.
+@test sim.grains[1].force[1] ≈ 0.
+@test sim.grains[1].force[2] ≈ 0.
+
+info("Wall present but no contact")
+sim = Granular.createSimulation()
+sim.ocean = Granular.createRegularOceanGrid([1, 1, 1], [10., 20., 1.0])
+Granular.addGrainCylindrical!(sim, [ 0., 0.], 1., 2., verbose=false)
+Granular.addWallLinearFrictionless!(sim, [1., 0.], +2.01, verbose=false)
+sim.walls[1].contact_viscosity_normal = 1e3
+Granular.setTimeStep!(sim, verbose=false)
+Granular.interactWalls!(sim)
+@test sim.walls[1].force ≈ 0.
+@test sim.grains[1].force[1] ≈ 0.
+@test sim.grains[1].force[2] ≈ 0.
+
+info("Wall at -x")
+sim = Granular.createSimulation()
+sim.ocean = Granular.createRegularOceanGrid([1, 1, 1], [10., 20., 1.0])
+Granular.addGrainCylindrical!(sim, [ 0., 0.], 1., 2., verbose=false)
+Granular.addWallLinearFrictionless!(sim, [1., 0.], -1. + .01, verbose=false)
+sim.walls[1].contact_viscosity_normal = 1e3
+Granular.setTimeStep!(sim, verbose=false)
+Granular.interactWalls!(sim)
+@test sim.walls[1].force < 0.
+@test sim.grains[1].force[1] > 0.
+@test sim.grains[1].force[2] ≈ 0.
+
+info("Wall at +x")
+sim = Granular.createSimulation()
+sim.ocean = Granular.createRegularOceanGrid([1, 1, 1], [10., 20., 1.0])
+Granular.addGrainCylindrical!(sim, [ 0., 0.], 1., 2., verbose=false)
+Granular.addWallLinearFrictionless!(sim, [1., 0.], +1. - .01, verbose=false)
+sim.walls[1].contact_viscosity_normal = 1e3
+Granular.setTimeStep!(sim, verbose=false)
+Granular.interactWalls!(sim)
+@test sim.walls[1].force > 0.
+@test sim.grains[1].force[1] < 0.
+@test sim.grains[1].force[2] ≈ 0.
+
+info("Wall at -y")
+sim = Granular.createSimulation()
+sim.ocean = Granular.createRegularOceanGrid([1, 1, 1], [10., 20., 1.0])
+Granular.addGrainCylindrical!(sim, [ 0., 0.], 1., 2., verbose=false)
+Granular.addWallLinearFrictionless!(sim, [0., 1.], -1. + .01, verbose=false)
+sim.walls[1].contact_viscosity_normal = 1e3
+Granular.setTimeStep!(sim, verbose=false)
+Granular.interactWalls!(sim)
+@test sim.walls[1].force < 0.
+@test sim.grains[1].force[1] ≈ 0.
+@test sim.grains[1].force[2] > 0.
+
+info("Wall at +y")
+sim = Granular.createSimulation()
+sim.ocean = Granular.createRegularOceanGrid([1, 1, 1], [10., 20., 1.0])
+Granular.addGrainCylindrical!(sim, [ 0., 0.], 1., 2., verbose=false)
+Granular.addWallLinearFrictionless!(sim, [0., 1.], +1. - .01, verbose=false)
+sim.walls[1].contact_viscosity_normal = 1e3
+Granular.setTimeStep!(sim, verbose=false)
+Granular.interactWalls!(sim)
+@test sim.walls[1].force > 0.
+@test sim.grains[1].force[1] ≈ 0.
+@test sim.grains[1].force[2] < 0.
+
+info("Full collision with wall")
+sim = Granular.createSimulation()
+sim.ocean = Granular.createRegularOceanGrid([1, 1, 1], [10., 20., 1.0])
+Granular.addGrainCylindrical!(sim, [1.2, 0.], 1., 2., verbose=false)
+Granular.addWallLinearFrictionless!(sim, [1., 0.], 0., verbose=false)
+sim.walls[1].contact_viscosity_normal = 1e3
+sim.grains[1].lin_vel[1] = -0.2
+Granular.setTimeStep!(sim, verbose=false)
+Granular.setTotalTime!(sim, 5.)
+lin_vel0 = sim.grains[1].lin_vel
+Granular.run!(sim)
+@test sim.grains[1].lin_vel[1] < abs(lin_vel0[1])
+@test sim.grains[1].lin_vel[2] ≈ 0.
+lin_vel1 = sim.grains[1].lin_vel
+
+sim = Granular.createSimulation()
+sim.ocean = Granular.createRegularOceanGrid([1, 1, 1], [10., 20., 1.0])
+Granular.addGrainCylindrical!(sim, [1.2, 0.], 1., 2., verbose=false)
+Granular.addWallLinearFrictionless!(sim, [1., 0.], 0., verbose=false)
+sim.walls[1].contact_viscosity_normal = 1e4
+sim.grains[1].lin_vel[1] = -0.2
+Granular.setTimeStep!(sim, verbose=false)
+Granular.setTotalTime!(sim, 5.)
+lin_vel0 = sim.grains[1].lin_vel
+Granular.run!(sim)
+@test sim.grains[1].lin_vel[1] < abs(lin_vel0[1])
+@test sim.grains[1].lin_vel[1] < lin_vel1[1]
+@test sim.grains[1].lin_vel[2] ≈ 0.
+
info("# Testing wall dynamics")
