commit 76b8123dd2aaa785b6e46792988949ac0122cd9a
parent bfb62ecae69cc80217011850ff9fce1839871fbb
Author: Anders Damsgaard <andersd@riseup.net>
Date: Fri, 10 Nov 2017 08:22:06 -0600
add flags to allow x or y movement although the grain is fixed
Diffstat:
5 files changed, 88 insertions(+), 12 deletions(-)
diff --git a/examples/shear.jl b/examples/shear.jl
@@ -16,6 +16,7 @@ const nx = 10
const ny = 50
const r_min = 0.2
const r_max = 1.0
+const z_thickness =
Granular.regularPacking!(sim, [nx, ny], r_min, r_max)
# Create a grid for contact searching spanning the extent of the grains
@@ -83,6 +84,7 @@ for grain in grains
y_top = grain.lin_pos[2] + grain.contact_radius
end
end
+# TODO
Granular.addDynamicWall!(sim, normal=[0., -1.], pos=y_top,
boundary_condition="normal stress",
normal_stress=N)
@@ -100,7 +102,7 @@ end
const fixed_thickness = 2. * r_max
for grain in grains
if grain.lin_pos[2] <= fixed_thickness
- grain.fixed = true
+ grain.fixed = true # set x and y acceleration to zero
end
end
@@ -127,10 +129,10 @@ sim_cons = deepcopy(sim)
################################################################################
# Select a shear velocity for the consolidation [m/s]
-const v_shear = 0.1
+const vel_shear = 0.1
# Rename the simulation so we don't overwrite output from the previous step
-sim.id = "$(id_prefix)-shear-N$(N)Pa-v_shear$(v_shear)m-s"
+sim.id = "$(id_prefix)-shear-N$(N)Pa-vel_shear$(vel_shear)m-s"
# Set all linear and rotational velocities to zero
Granular.zeroKinematics!(sim)
@@ -146,7 +148,7 @@ for grain in grains
# overrides the `fixed` flag
grain.allow_y_acc = true
- grain.lin_vel[1] = v_shear
+ grain.lin_vel[1] = vel_shear
end
end
@@ -156,9 +158,54 @@ Granular.resetTime!(sim)
# Set the simulation time to run the shear experiment for
Granular.setTotalTime!(sim, 15.0)
-# Run the consolidation experiment
-Granular.run!(sim)
+# Run the shear experiment
+time = Float64[]
+shear_stress = Float64[]
+shear_strain = Float64[]
+dilation = Float64[]
+const thickness_initial = sim.walls[1].pos - y_bot
+x_min = +Inf
+x_max = -Inf
+for grain in grains
+ if x_min > grain.lin_pos[1] - grain.contact_radius
+ x_min = grain.lin_pos[1] - grain.contact_radius
+ end
+ if x_max < grain.lin_pos[1] + grain.contact_radius
+ x_max = grain.lin_pos[1] + grain.contact_radius
+ end
+end
+const surface_area = (x_max - x_min)
+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
+
+ append!(time, sim.time)
+
+ # Determine the current shear stress
+ for grain in sim.grains
+ if grain.fixed && grain.allow_y_acc
+ shear_force += -grain.force[1]
+ end
+ end
+ append!(shear_stress, shear_force/surface_area)
+
+ # Determine the current shear strain
+ append!(shear_strain, sim.time*vel_shear/thickness_initial)
+
+ # Determine the current dilation
+ append!(dilation, (sim.walls[1].pos - y_bot)/thickness_initial)
+
# Save the simulation state to disk in case we need to reuse the sheared state
Granular.writeSimulation(sim)
+# Plot time vs. shear stress and dilation
+# TODO
+
+# Plot shear strain vs. shear stress and dilation
+# TODO
+
+# Plot time vs. shear strain (boring when the shear experiment is rate controlled)
+# TODO
diff --git a/src/datatypes.jl b/src/datatypes.jl
@@ -30,6 +30,8 @@ mutable struct GrainCylindrical
# Kinematic constraint flags
fixed::Bool
+ allow_x_acc::Bool
+ allow_y_acc::Bool
rotating::Bool
enabled::Bool
@@ -99,6 +101,8 @@ mutable struct GrainArrays
# Kinematic constraint flags
fixed::Vector{Int}
+ allow_x_acc::Vector{Int}
+ allow_y_acc::Vector{Int}
rotating::Vector{Int}
enabled::Vector{Int}
diff --git a/src/grain.jl b/src/grain.jl
@@ -25,7 +25,9 @@ export addGrainCylindrical!
ocean_drag_coeff_horiz,
atmosphere_drag_coeff_vert,
atmosphere_drag_coeff_horiz,
- pressure, fixed, rotating, enabled, verbose,
+ pressure, fixed,
+ allow_x_acc, allow_y_acc,
+ rotating, enabled, verbose,
ocean_grid_pos, atmosphere_grid_pos,
n_contact, granular_stress, ocean_stress,
atmosphere_stress])
@@ -84,6 +86,8 @@ are optional, and come with default values. The only required arguments are
for atmosphere against grain bottom [-].
* `pressure::Float64 = 0.`: current compressive stress on grain [Pa].
* `fixed::Bool = false`: grain is fixed to a constant velocity (e.g. zero).
+* `allow_x_acc::Bool = false`: override `fixed` along `x`.
+* `allow_y_acc::Bool = false`: override `fixed` along `y`.
* `rotating::Bool = true`: grain is allowed to rotate.
* `enabled::Bool = true`: grain interacts with other grains.
* `verbose::Bool = true`: display diagnostic information during the function
@@ -162,6 +166,8 @@ function addGrainCylindrical!(simulation::Simulation,
# H&C2011
pressure::Float64 = 0.,
fixed::Bool = false,
+ allow_x_acc::Bool = false,
+ allow_y_acc::Bool = false,
rotating::Bool = true,
enabled::Bool = true,
verbose::Bool = true,
@@ -230,6 +236,8 @@ function addGrainCylindrical!(simulation::Simulation,
torque,
fixed,
+ allow_x_acc,
+ allow_y_acc,
rotating,
enabled,
@@ -347,6 +355,8 @@ function convertGrainDataToArrays(simulation::Simulation)
Array{Int}(length(simulation.grains)),
Array{Int}(length(simulation.grains)),
Array{Int}(length(simulation.grains)),
+ Array{Int}(length(simulation.grains)),
+ Array{Int}(length(simulation.grains)),
Array{Float64}(length(simulation.grains)),
Array{Float64}(length(simulation.grains)),
@@ -399,6 +409,8 @@ function convertGrainDataToArrays(simulation::Simulation)
ifarr.torque[3, i] = simulation.grains[i].torque
ifarr.fixed[i] = Int(simulation.grains[i].fixed)
+ ifarr.allow_x_acc[i] = Int(simulation.grains[i].allow_x_acc)
+ ifarr.allow_y_acc[i] = Int(simulation.grains[i].allow_y_acc)
ifarr.rotating[i] = Int(simulation.grains[i].rotating)
ifarr.enabled[i] = Int(simulation.grains[i].enabled)
@@ -470,6 +482,8 @@ function deleteGrainArrays!(ifarr::GrainArrays)
ifarr.torque = f2
ifarr.fixed = i1
+ ifarr.allow_x_acc = i1
+ ifarr.allow_y_acc = i1
ifarr.rotating = i1
ifarr.enabled = i1
@@ -529,9 +543,11 @@ function printGrainInfo(f::GrainCylindrical)
println(" ang_acc: $(f.ang_acc) rad/s^2")
println(" torque: $(f.torque) N*m\n")
- println(" fixed: $(f.fixed)")
- println(" rotating: $(f.rotating)")
- println(" enabled: $(f.enabled)\n")
+ println(" fixed: $(f.fixed)")
+ println(" allow_x_acc: $(f.allow_x_acc)")
+ println(" allow_y_acc: $(f.allow_y_acc)")
+ println(" rotating: $(f.rotating)")
+ println(" enabled: $(f.enabled)\n")
println(" k_n: $(f.contact_stiffness_normal) N/m")
println(" k_t: $(f.contact_stiffness_tangential) N/m")
diff --git a/src/temporal.jl b/src/temporal.jl
@@ -171,7 +171,10 @@ export resetTime!
resetTime!(simulation)
Reset the current time to zero, and reset output file counters in order to
-restart a simulation.
+restart a simulation. This function does not overwrite the time step
+(`Simulation.time_step`), the output
+file interval (`Simulation.file_time_step`), or the total simulation time
+(`Simulation.time_total`).
# Arguments
* `simulation::Simulation`: the simulation object for which to reset the
@@ -183,3 +186,4 @@ function resetTime!(sim::Simulation)
sim.file_number = 0
sim.file_time_since_output_file = 0.
end
+
diff --git a/src/temporal_integration.jl b/src/temporal_integration.jl
@@ -54,7 +54,12 @@ function updateGrainKinematicsTwoTermTaylor!(grain::GrainCylindrical,
grain.ang_acc = grain.torque/grain.moment_of_inertia
if grain.fixed
- fill!(grain.lin_acc, 0.)
+ if !grain.allow_x_acc
+ grain.lin_acc[1] = 0.
+ end
+ if !grain.allow_y_acc
+ grain.lin_acc[2] = 0.
+ end
grain.ang_acc = 0.
elseif !grain.rotating
grain.ang_acc = 0.
