Granular.jl

commit 6e8adb0d8a14aa77aaa516f197c72bbe242a9fbb
parent ccd9921a7fd354df82bcfcfb115844cc26278beb
Author: Anders Damsgaard <andersd@riseup.net>
Date:   Tue, 14 Nov 2017 16:12:33 -0500

first implementaiton of grain-wall interaction

Diffstat:
M
src/interaction.jl
 | 
42
++++++++++++++++++++++++++++++++++++++++++
M
src/simulation.jl
 | 
5
++++-
M
src/wall.jl
 | 
15
+++++++++------

3 files changed, 55 insertions(+), 7 deletions(-)
diff --git a/src/interaction.jl b/src/interaction.jl
@@ -5,6 +5,9 @@ export interact!
     interact!(simulation::Simulation)
 
 Resolve mechanical interaction between all particle pairs.
+
+# Arguments
+* `simulation::Simulation`: the simulation object containing the grains.
 """
 function interact!(simulation::Simulation)
     for i=1:length(simulation.grains)
@@ -28,6 +31,45 @@ function interact!(simulation::Simulation)
     nothing
 end
 
+
+"""
+    interactWalls!(sim)
+
+Find and resolve interactions between the dynamic walls (`simulation.walls`) and
+the grains.  The contact model uses linear elasticity, with stiffness based on
+the grain Young's modulus `grian.E` or elastic stiffness `grain.k_n`.  The
+interaction is frictionless in the tangential direction.
+
+# Arguments
+* `simulation::Simulation`: the simulation object containing the grains and
+    dynamic walls.
+"""
+function interactWalls!(sim::Simulation)
+
+    δ_n::Float64 = 0.0
+    k_n::Float64 = 0.0
+
+    for iw=1:length(sim.walls)
+        for i=1:length(sim.grains)
+
+            # get overlap distance by projecting grain position onto wall-normal
+            # vector
+            δ_n = dot(sim.walls[iw].normal, sim.grains[i].lin_pos) -
+                sim.walls[iw].pos - sim.grains[i].contact_radius
+
+            if sim.grains[i].youngs_modulus > 0.
+                k_n = sim.grains[i].youngs_modulus *
+                    sim.grains[i].thickness
+            else
+                k_n = sim.grains[i].contact_stiffness_normal
+            end
+
+            sim.walls[iw].force += k_n * abs(δ_n)
+            sim.grains[i].force += k_n * abs(δ_n) * sim.grains[iw].normal
+        end
+    end
+end
+
 export interactGrains!
 """
     interactGrains!(simulation::Simulation, i::Int, j::Int, ic::Int)
diff --git a/src/simulation.jl b/src/simulation.jl
@@ -248,13 +248,16 @@ function disableGrain!(simulation::Simulation, i::Int)
 end
 
 export zeroForcesAndTorques!
-"Sets the `force` and `torque` values of all grains to zero."
+"Sets the `force` and `torque` values of all grains and dynamic walls to zero."
 function zeroForcesAndTorques!(simulation::Simulation)
     for grain in simulation.grains
         grain.force .= grain.external_body_force
         grain.torque = 0.
         grain.pressure = 0.
     end
+    for wall in simulation.walls
+        wall.force = 0.
+    end
     nothing
 end
 
diff --git a/src/wall.jl b/src/wall.jl
@@ -13,7 +13,9 @@ The only required arguments are
 # Arguments
 * `simulation::Simulation`: the simulation object where the wall should be
     added to.
-* `normal::Vector{Float64}`: 2d vector denoting the normal to the wall [m].
+* `normal::Vector{Float64}`: 2d vector denoting the normal to the wall [m].  The
+    wall will only interact in the opposite direction of this vector, so the
+    normal vector should point in the direction of the grains.
 * `pos::Float64`: position along axis parallel to the normal vector [m].
 * `bc::String="fixed"`: boundary condition, possible values are `"fixed"`
     (default), `"normal stress"`, or `"velocity"`.
@@ -52,11 +54,11 @@ Granular.addWallLinearFrictionless!(sim, [0., 1.], 1.5,
                                     vel=-0.5)
 ```
 
-To create a wall parallel to the *x* axis with a constant normal stress of 100
-kPa:
+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, [1., 0.], 3.5,
+Granular.addWallLinearFrictionless!(sim, [0., -1.], 3.5,
                                     bc="normal stress",
                                     normal_stress=100e3)
 ```
@@ -78,11 +80,12 @@ function addWallLinearFrictionless!(simulation::Simulation,
               "$normal)")
     end
 
-    if !(normal ≈ [1., 0.]) && !(normal ≈ [0., 1.])
+    if !(normal ≈ [1., 0.]) && !(normal ≈ [0., 1.]) &&
+        !(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.