commit 31b70e6d658735acfb61396662856f895c8354cb
parent 8d576c42b7df677e8084e01fa5be310bad08148c
Author: Anders Damsgaard <andersd@riseup.net>
Date: Thu, 16 Nov 2017 12:18:31 -0800
fix wall initialization, update grain size in example
Diffstat:
3 files changed, 31 insertions(+), 28 deletions(-)
diff --git a/examples/shear.jl b/examples/shear.jl
@@ -19,8 +19,8 @@ const nx = 10 # Grains along x (horizontal)
const ny = 50 # Grains along y (vertical)
# Grain-size parameters
-const r_min = 0.2 # Min. grain radius [m]
-const r_max = 1.0 # Max. grain radius [m]
+const r_min = 0.01 # Min. grain radius [m]
+const r_max = 0.1 # Max. grain radius [m]
const gsd_type = "powerlaw" # "powerlaw" or "uniform" sizes between r_min and r_max
const gsd_powerlaw_exponent = -1.8 # GSD power-law exponent
const gsd_seed = 1 # Value to seed random-size generation
@@ -36,7 +36,7 @@ const rotating = true # Allow grain rotation
const N = 10e3
# Shear velocity to apply to the top grains [m/s]
-const vel_shear = 0.1
+const vel_shear = 0.01
################################################################################
#### Step 1: Create a loose granular assemblage and let it settle at -y #
@@ -66,10 +66,13 @@ Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "north south",
verbose=false)
Granular.setGridBoundaryConditions!(sim.ocean, "periodic", "east west")
-# Add gravitational acceleration to all grains and disable ocean-grid drag
+# Add gravitational acceleration to all grains and disable ocean-grid drag.
+# Also add viscous energy dissipation between grains, which is disabled before
+# consolidation and shear.
for grain in sim.grains
Granular.addBodyForce!(grain, grain.mass*g)
Granular.disableOceanDrag!(grain)
+ grain.contact_viscosity_normal = 1e4 # N/(m/s)
end
# Automatically set the computational time step based on grain sizes and
@@ -77,13 +80,14 @@ end
Granular.setTimeStep!(sim)
# Set the total simulation time for this step [s]
-Granular.setTotalTime!(sim, 30.)
+# This value may need tweaking if grain sizes or numbers are adjusted.
+Granular.setTotalTime!(sim, 1.)
# Set the interval in model time between simulation files [s]
-Granular.setOutputFileInterval!(sim, .2)
+Granular.setOutputFileInterval!(sim, .01)
# Visualize the grain-size distribution
-#Granular.plotGrainSizeDistribution(sim)
+Granular.plotGrainSizeDistribution(sim)
# Start the simulation
Granular.run!(sim)
@@ -112,9 +116,10 @@ Granular.zeroKinematics!(sim)
# Add a dynamic wall to the top which adds a normal stress downwards. The
# normal of this wall is downwards, and we place it at the top of the granular
-# assemblage
+# assemblage. Here, the inter-grain viscosity is also removed.
y_top = -Inf
for grain in sim.grains
+ grain.contact_viscosity_normal = 0.
if y_top < grain.lin_pos[2] + grain.contact_radius
y_top = grain.lin_pos[2] + grain.contact_radius
end
@@ -144,7 +149,7 @@ end
Granular.resetTime!(sim)
# Set the simulation time to run the consolidation for
-Granular.setTotalTime!(sim, 5.0)
+Granular.setTotalTime!(sim, 1.0)
# Run the consolidation experiment, and monitor top wall position over time
time = Float64[]
@@ -155,6 +160,8 @@ while sim.time < sim.time_total
for i=1:100 # run for 100 steps before measuring shear stress and dilation
Granular.run!(sim, single_step=true)
end
+ println(sim.walls[1])
+ println(sim.walls[1].pos)
append!(time, sim.time)
append!(compaction, sim.walls[1].pos)
@@ -203,7 +210,7 @@ Granular.zeroKinematics!(sim)
Granular.resetTime!(sim)
# Set the simulation time to run the shear experiment for
-Granular.setTotalTime!(sim, 15.0)
+Granular.setTotalTime!(sim, 2.0)
# Run the shear experiment
time = Float64[]
diff --git a/src/temporal_integration.jl b/src/temporal_integration.jl
@@ -177,14 +177,12 @@ function updateWallKinematicsTwoTermTaylor!(wall::WallLinearFrictionless,
simulation::Simulation)
if wall.bc == "fixed"
return nothing
- end
-
- if wall.bc == "velocity"
+ elseif wall.bc == "velocity"
wall.acc = 0.0
+ elseif wall.bc == "normal stress"
+ wall.acc = (wall.force + wall.normal_stress*wall.surface_area)/wall.mass
else
- # Normal force directed along normal
- f_n::Float64 = -wall.normal_stress*wall.surface_area
- wall.acc = (wall.force + f_n)/wall.mass
+ error("wall boundary condition was not understood ($(wall.bc))")
end
wall.pos +=
@@ -213,14 +211,12 @@ function updateWallKinematicsThreeTermTaylor!(wall::WallLinearFrictionless,
if wall.bc == "fixed"
return nothing
- end
-
- # Normal load = normal stress times wall surface area, directed along normal
-
- if wall.bc == "velocity"
+ elseif wall.bc == "velocity"
wall.acc = 0.0
- else
+ elseif wall.bc == "normal stress"
wall.acc = (wall.force + wall.normal_stress*wall.surface_area)/wall.mass
+ else
+ error("wall boundary condition was not understood ($(wall.bc))")
end
# Temporal gradient in acceleration using backwards differences
diff --git a/src/wall.jl b/src/wall.jl
@@ -67,9 +67,9 @@ function addWallLinearFrictionless!(simulation::Simulation,
normal::Vector{Float64},
pos::Float64;
bc::String = "fixed",
- mass::Float64 = NaN,
- thickness::Float64 = NaN,
- surface_area::Float64 = NaN,
+ mass::Float64 = -1.,
+ thickness::Float64 = -1.,
+ surface_area::Float64 = -1.,
normal_stress::Float64 = 0.,
vel::Float64 = 0.,
acc::Float64 = 0.,
@@ -94,7 +94,7 @@ function addWallLinearFrictionless!(simulation::Simulation,
end
# if not set, set wall mass to equal the mass of all grains.
- if isnan(mass)
+ if mass < 0.
if length(simulation.grains) < 1
error("If wall mass is not specified, walls should be added " *
"after grains have been added to the simulation.")
@@ -109,7 +109,7 @@ function addWallLinearFrictionless!(simulation::Simulation,
end
# if not set, set wall thickness to equal largest grain thickness
- if isnan(thickness)
+ if thickness < 0.
if length(simulation.grains) < 1
error("If wall thickness is not specified, walls should be added " *
"after grains have been added to the simulation.")
@@ -126,7 +126,7 @@ function addWallLinearFrictionless!(simulation::Simulation,
end
# if not set, set wall surface area from the ocean grid
- if isnan(surface_area) && bc != "fixed"
+ if surface_area < 0. && bc != "fixed"
if typeof(simulation.ocean.input_file) == Bool
error("simulation.ocean must be set beforehand")
end
