Granular.jl

Julia package for granular dynamics simulation
git clone git://src.adamsgaard.dk/Granular.jl
Log | Files | Refs | README | LICENSE

commit 5951464686ab69926cfb9d4905f42ea898d38789
parent f2a496feac2b8cbb587385d9a4068a508a4f832e
Author: Anders Damsgaard <andersd@riseup.net>
Date:   Wed, 15 Nov 2017 16:02:34 -0500

wall normal must be positive, add wall temporal integration, add more tests

Diffstat:
Mexamples/shear.jl | 1-
Msrc/datatypes.jl | 2++
Msrc/interaction.jl | 7+------
Msrc/temporal_integration.jl | 120+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--
Msrc/wall.jl | 37++++++++++++++++++++++++++++++++-----
Mtest/wall.jl | 223++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-----
6 files changed, 364 insertions(+), 26 deletions(-)

diff --git a/examples/shear.jl b/examples/shear.jl @@ -16,7 +16,6 @@ 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 diff --git a/src/datatypes.jl b/src/datatypes.jl @@ -77,8 +77,10 @@ mutable struct WallLinearFrictionless bc::String # Boundary condition mass::Float64 # Mass, used when bc != "fixed" thickness::Float64 # Wall thickness + surface_area::Float64 # Wall surface area normal_stress::Float64 # Normal stress when bc == "normal stress" vel::Float64 # Velocity (constant when bc == "normal stress") + acc::Float64 # Acceleration (zero when bc == "velocity") force::Float64 # Sum of normal forces on wall end diff --git a/src/interaction.jl b/src/interaction.jl @@ -57,11 +57,6 @@ function interactWalls!(sim::Simulation) δ_n = abs(dot(sim.walls[iw].normal, sim.grains[i].lin_pos) - sim.walls[iw].pos) - sim.grains[i].contact_radius - #println(dot(sim.walls[iw].normal, sim.grains[i].lin_pos)) - #println(sim.walls[iw].pos) - #println(sim.grains[i].contact_radius) - #println(δ_n) - if δ_n < 0. if sim.grains[i].youngs_modulus > 0. k_n = sim.grains[i].youngs_modulus * sim.grains[i].thickness @@ -69,8 +64,8 @@ function interactWalls!(sim::Simulation) k_n = sim.grains[i].contact_stiffness_normal end - sim.walls[iw].force += -k_n * δ_n sim.grains[i].force += k_n * abs(δ_n) * sim.walls[iw].normal + sim.walls[iw].force += k_n * δ_n end end end diff --git a/src/temporal_integration.jl b/src/temporal_integration.jl @@ -46,7 +46,7 @@ end export updateGrainKinematicsTwoTermTaylor! """ Use a two-term Taylor expansion for integrating the kinematic degrees of freedom -for an `grain`. +for a `grain`. """ function updateGrainKinematicsTwoTermTaylor!(grain::GrainCylindrical, simulation::Simulation) @@ -80,7 +80,7 @@ end export updateGrainKinematicsThreeTermTaylor! """ Use a three-term Taylor expansion for integrating the kinematic degrees of -freedom for an `grain`. +freedom for a `grain`. """ function updateGrainKinematicsThreeTermTaylor!(grain::GrainCylindrical, simulation::Simulation) @@ -122,3 +122,119 @@ function updateGrainKinematicsThreeTermTaylor!(grain::GrainCylindrical, 0.5 * d_ang_acc_dt * simulation.time_step^2. nothing end + +export updateWallKinematics! +""" + updateWallKinematics!(simulation::Simulation[, + method::String = "Three-term Taylor"]) + +Update the wall kinematic parameters using a temporal integration scheme, +the current force and torque balance, and gravitational acceleration. If the +simulation contains a grid with periodic boundaries, affected wall positions +are adjusted accordingly. + +# Arguments +* `simulation::Simulation`: update the wall positions in this object + according to temporal integration of length `simulation.time_step`. +* `method::String = "Three-term Taylor"`: the integration method to use. + Available methods are "Two-term Taylor" and "Three-term Taylor". The + three-term Taylor expansion is slightly more computationally expensive than + the two-term Taylor expansion, but offers an order-of-magnitude increase in + precision of wall positions. The two-term expansion can obtain similar + precision if the time step is 1/10 the length. +""" +function updateWallKinematics!(simulation::Simulation; + method::String = "Three-term Taylor") + + if method == "Two-term Taylor" + for wall in simulation.walls + if wall.bc == "fixed" + continue + end + updateWallKinematicsTwoTermTaylor!(wall, simulation) + end + elseif method == "Three-term Taylor" + for wall in simulation.walls + if wall.bc == "fixed" + continue + end + updateWallKinematicsThreeTermTaylor!(wall, simulation) + end + else + error("Unknown integration method '$method'") + end + nothing +end + +export updateWallKinematicsTwoTermTaylor! +""" + updateWallKinematicsTwoTermTaylor!(wall, simulation) + +Use a two-term Taylor expansion for integrating the kinematic degrees of freedom +for a `wall`. +""" +function updateWallKinematicsTwoTermTaylor!(wall::WallLinearFrictionless, + simulation::Simulation) + if wall.bc == "fixed" + return nothing + end + + if wall.bc == "velocity" + wall.acc = 0.0 + else + # Normal force directed along normal + f_n::Float64 = -wall.normal_stress*wall.surface_area + wall.acc = (wall.force + f_n)/wall.mass + end + + wall.pos += + wall.vel * simulation.time_step + + 0.5*wall.acc * simulation.time_step^2.0 + + wall.vel += wall.acc * simulation.time_step + nothing +end + + +export updateWallKinematicsThreeTermTaylor! +""" + updateWallKinematicsThreeTermTaylor!(wall, simulation) + +Use a two-term Taylor expansion for integrating the kinematic degrees of freedom +for a `wall`. +""" +function updateWallKinematicsThreeTermTaylor!(wall::WallLinearFrictionless, + simulation::Simulation) + if simulation.time_iteration == 0 + acc_0 = 0. + else + acc_0 = wall.acc + end + + if wall.bc == "fixed" + return nothing + end + + # Normal load = normal stress times wall surface area, directed along normal + + if wall.bc == "velocity" + wall.acc = 0.0 + else + f_n::Float64 = -wall.normal_stress*wall.surface_area + wall.acc = (wall.force + f_n)/wall.mass + end + + # Temporal gradient in acceleration using backwards differences + d_acc_dt = (wall.acc - acc_0)/simulation.time_step + + wall.pos += + wall.vel * simulation.time_step + + 0.5*wall.acc * simulation.time_step^2.0 + + 1.0/6.0 * d_acc_dt * simulation.time_step^3.0 + + wall.vel += wall.acc * simulation.time_step + + 0.5 * d_acc_dt * simulation.time_step^2.0 + + nothing +end + diff --git a/src/wall.jl b/src/wall.jl @@ -3,7 +3,7 @@ export addWallLinearFrictionless! """ function addWallLinear!(simulation, normal, pos[, bc, mass, thickness, - normal_stress, vel, force, verbose]) + normal_stress, vel, acc, force, verbose]) Creates and adds a linear (flat) and frictionless dynamic wall to a grain to a simulation. Most of the arguments are optional, and come with default values. @@ -58,7 +58,7 @@ To create a wall parallel to the *y* axis pushing downwards with a constant normal stress of 100 kPa, starting at a position of y = 3.5 m: ```julia -Granular.addWallLinearFrictionless!(sim, [0., -1.], 3.5, +Granular.addWallLinearFrictionless!(sim, [0., 1.], 3.5, bc="normal stress", normal_stress=100e3) ``` @@ -69,8 +69,10 @@ function addWallLinearFrictionless!(simulation::Simulation, bc::String = "fixed", mass::Float64 = NaN, thickness::Float64 = NaN, + surface_area::Float64 = NaN, normal_stress::Float64 = 0., vel::Float64 = 0., + acc::Float64 = 0., force::Float64 = 0., verbose::Bool=true) @@ -80,12 +82,15 @@ function addWallLinearFrictionless!(simulation::Simulation, "$normal)") end - if !(normal ≈ [1., 0.]) && !(normal ≈ [0., 1.]) && - !(normal ≈ [-1., 0.]) && !(normal ≈ [0., -1.]) + if bc != "fixed" && bc != "velocity" && bc != "normal stress" + error("Wall BC must be 'fixed', 'velocity', or 'normal stress'.") + end + + if !(normal ≈ [1., 0.]) && !(normal ≈ [0., 1.]) error("Currently only walls with normals orthogonal to the " * "coordinate system are allowed, i.e. normals parallel to the " * "x or y axes. Accepted values for `normal` " * - "are [±1., 0.] and [0., ±1.]. The passed normal was $normal") + "are [1., 0.] and [0., 1.]. The passed normal was $normal") end # if not set, set wall mass to equal the mass of all grains. @@ -120,14 +125,24 @@ function addWallLinearFrictionless!(simulation::Simulation, end end + # if not set, set wall surface area from the ocean grid + if isnan(surface_area) && bc != "fixed" + if typeof(simulation.ocean.input_file) == Bool + error("simulation.ocean must be set beforehand") + end + surface_area = getWallSurfaceArea(simulation, normal, thickness) + end + # Create wall object wall = WallLinearFrictionless(normal, pos, bc, mass, thickness, + surface_area, normal_stress, vel, + acc, force) # Add to simulation object @@ -158,6 +173,18 @@ function getWallSurfaceArea(sim::Simulation, wall_index::Integer) end nothing end +function getWallSurfaceArea(sim::Simulation, normal::Vector{Float64}, + thickness::Float64) + + if normal ≈ [1., 0.] + return (sim.ocean.yq[end,end] - sim.ocean.yq[1,1]) * thickness + elseif normal ≈ [0., 1.] + return (sim.ocean.xq[end,end] - sim.ocean.xq[1,1]) * thickness + else + error("Wall normal not understood") + end + nothing +end export getWallNormalStress """ diff --git a/test/wall.jl b/test/wall.jl @@ -4,6 +4,7 @@ info("#### $(basename(@__FILE__)) ####") +info("# Test wall initialization") info("Testing argument value checks") sim = Granular.createSimulation() Granular.addGrainCylindrical!(sim, [ 0., 0.], 10., 2., verbose=false) @@ -16,6 +17,9 @@ Granular.addGrainCylindrical!(sim, [ 0., 0.], 10., 2., verbose=false) @test_throws ErrorException Granular.addWallLinearFrictionless!(sim, [.1, .1, .1], 1.) +@test_throws ErrorException Granular.addWallLinearFrictionless!(sim, + [0., 1.], 1., + bc="asdf") sim = Granular.createSimulation() @test_throws ErrorException Granular.addWallLinearFrictionless!(sim, [1., 0.], 1.) @@ -46,64 +50,259 @@ 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") -info("...wall present but no contact") +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) +Granular.setTimeStep!(sim) 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") +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) +Granular.addWallLinearFrictionless!(sim, [1., 0.], +2.01, verbose=false) +Granular.setTimeStep!(sim) 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") +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) +Granular.setTimeStep!(sim) Granular.interactWalls!(sim) @test sim.walls[1].force > 0. -@test sim.grains[1].force[1] > 0. +@test sim.grains[1].force[1] < 0. @test sim.grains[1].force[2] ≈ 0. -info("...wall at +x") +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) +Granular.addWallLinearFrictionless!(sim, [1., 0.], +1. - .01, verbose=false) +Granular.setTimeStep!(sim) 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") +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) +Granular.setTimeStep!(sim) Granular.interactWalls!(sim) -@test sim.walls[1].force > 0. +@test sim.walls[1].force < 0. @test sim.grains[1].force[1] ≈ 0. @test sim.grains[1].force[2] > 0. -info("...wall at +y") +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) +Granular.addWallLinearFrictionless!(sim, [0., 1.], +1. - .01, verbose=false) +Granular.setTimeStep!(sim) Granular.interactWalls!(sim) @test sim.walls[1].force > 0. @test sim.grains[1].force[1] ≈ 0. @test sim.grains[1].force[2] < 0. + + +info("# Testing wall dynamics") + +info("Wall present, no contact, fixed (default)") +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) +Granular.setTimeStep!(sim) +Granular.interactWalls!(sim) +Granular.updateWallKinematics!(sim) +@test sim.walls[1].force ≈ 0. +@test sim.walls[1].acc ≈ 0. +@test sim.walls[1].vel ≈ 0. +@test sim.walls[1].pos ≈ -1.01 +@test sim.grains[1].force[1] ≈ 0. +@test sim.grains[1].force[2] ≈ 0. + +info("Wall present, no contact, fixed (TY2)") +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) +Granular.setTimeStep!(sim) +Granular.interactWalls!(sim) +Granular.updateWallKinematics!(sim, method="Two-term Taylor") +@test sim.walls[1].force ≈ 0. +@test sim.walls[1].acc ≈ 0. +@test sim.walls[1].vel ≈ 0. +@test sim.walls[1].pos ≈ -1.01 +@test sim.grains[1].force[1] ≈ 0. +@test sim.grains[1].force[2] ≈ 0. + +info("Wall present, no contact, fixed (TY3)") +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) +Granular.setTimeStep!(sim) +Granular.interactWalls!(sim) +@test_throws ErrorException Granular.updateWallKinematics!(sim, method="asdf") +Granular.updateWallKinematics!(sim, method="Three-term Taylor") +@test sim.walls[1].force ≈ 0. +@test sim.walls[1].acc ≈ 0. +@test sim.walls[1].vel ≈ 0. +@test sim.walls[1].pos ≈ -1.01 +@test sim.grains[1].force[1] ≈ 0. +@test sim.grains[1].force[2] ≈ 0. + +info("Wall present, contact, fixed") +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) +Granular.setTimeStep!(sim) +Granular.interactWalls!(sim) +Granular.updateWallKinematics!(sim) +@test sim.walls[1].acc ≈ 0. +@test sim.walls[1].vel ≈ 0. +@test sim.walls[1].pos ≈ -1.01 + +info("Wall present, no contact, velocity BC") +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, + bc="velocity", vel=1.0, + verbose=false) +Granular.setTimeStep!(sim) +Granular.interactWalls!(sim) +Granular.updateWallKinematics!(sim) +@test sim.walls[1].force ≈ 0. +@test sim.walls[1].acc ≈ 0. +@test sim.walls[1].vel ≈ 1. +@test sim.walls[1].pos > -1.01 +@test sim.grains[1].force[1] ≈ 0. +@test sim.grains[1].force[2] ≈ 0. + +info("Wall present, no contact, velocity BC (TY2)") +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, + bc="velocity", vel=1.0, + verbose=false) +Granular.setTimeStep!(sim) +Granular.interactWalls!(sim) +Granular.updateWallKinematics!(sim, method="Two-term Taylor") +@test sim.walls[1].force ≈ 0. +@test sim.walls[1].acc ≈ 0. +@test sim.walls[1].vel ≈ 1. +@test sim.walls[1].pos > -1.01 +@test sim.grains[1].force[1] ≈ 0. +@test sim.grains[1].force[2] ≈ 0. + +info("Wall present, no contact, velocity BC (TY3)") +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, + bc="velocity", vel=1.0, + verbose=false) +Granular.setTimeStep!(sim) +Granular.interactWalls!(sim) +@test_throws ErrorException Granular.updateWallKinematics!(sim, method="asdf") +Granular.updateWallKinematics!(sim, method="Three-term Taylor") +@test sim.walls[1].force ≈ 0. +@test sim.walls[1].acc ≈ 0. +@test sim.walls[1].vel ≈ 1. +@test sim.walls[1].pos > -1.01 +@test sim.grains[1].force[1] ≈ 0. +@test sim.grains[1].force[2] ≈ 0. + +info("Wall present, contact, velocity BC (TY2)") +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.], -0.9, + bc="velocity", vel=1.0, + verbose=false) +Granular.setTimeStep!(sim) +Granular.interactWalls!(sim) +Granular.updateWallKinematics!(sim, method="Two-term Taylor") +@test sim.walls[1].bc == "velocity" +@test sim.walls[1].acc ≈ 0. +@test sim.walls[1].vel ≈ 1. +@test sim.walls[1].pos > -0.9 + +info("Wall present, contact, velocity BC (TY2)") +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.], -0.9, + bc="velocity", vel=1.0, + verbose=false) +Granular.setTimeStep!(sim) +Granular.interactWalls!(sim) +Granular.updateWallKinematics!(sim, method="Two-term Taylor") +@test sim.walls[1].acc ≈ 0. +@test sim.walls[1].vel ≈ 1. +@test sim.walls[1].pos > -0.9 + +info("Wall present, contact, normal stress BC") +sim = Granular.createSimulation() +sim.ocean = Granular.createRegularOceanGrid([1, 1, 1], [2., 2., 1.]) +Granular.addGrainCylindrical!(sim, [ 1., 1.], 1., 1., verbose=false) +Granular.addWallLinearFrictionless!(sim, [1., 0.], 2., + bc="normal stress", + normal_stress=0., + verbose=false) +Granular.setTimeStep!(sim) +Granular.interactWalls!(sim) +Granular.updateWallKinematics!(sim) +@test sim.walls[1].force ≈ 0. +@test sim.walls[1].acc ≈ 0. +@test sim.walls[1].vel ≈ 0. +@test sim.grains[1].force[1] ≈ 0. +@test sim.grains[1].force[2] ≈ 0. + +info("Wall present, contact, normal stress BC") +sim = Granular.createSimulation() +sim.ocean = Granular.createRegularOceanGrid([1, 1, 1], [2., 2., 1.]) +Granular.addGrainCylindrical!(sim, [ 1., 1.], 1., 1., verbose=false) +sim.grains[1].fixed = true +Granular.addWallLinearFrictionless!(sim, [1., 0.], 2., + bc="normal stress", + normal_stress=1e4, + verbose=false) +Granular.setTimeStep!(sim) +Granular.interactWalls!(sim) +Granular.updateWallKinematics!(sim) +@test sim.walls[1].force ≈ 0. +@test sim.walls[1].acc < 0. +@test sim.walls[1].vel < 0. +@test sim.grains[1].force[1] ≈ 0. +@test sim.grains[1].force[2] ≈ 0. +for i=1:5 + Granular.interactWalls!(sim) + Granular.updateWallKinematics!(sim) + println(sim.walls[1].pos) +end +@test sim.walls[1].force < 0. +@test sim.walls[1].acc < 0. +@test sim.walls[1].vel < 0. +@test sim.walls[1].pos < 1. +@test sim.grains[1].force[1] > 0. +@test sim.grains[1].force[2] ≈ 0. +