timplement bond-parallel interaction - slidergrid - grid of elastic sliders on a frictional surface

git clone git://src.adamsgaard.dk/slidergrid

Log

Files

Refs

README

LICENSE

---

commit de268460b46629fdcfcfd155a1699fd3bd69b5c0

parent 73fcb298a3755f169f631813de24a4e2904510d3

Author: Anders Damsgaard 

Date:   Wed, 16 Mar 2016 11:41:29 -0700

implement bond-parallel interaction

Diffstat:
  M slider.c                            |      78 +++++++++++++++++++++++++------
  M slider.h                            |       5 +++--
  M vector_math.h                       |       2 +-

3 files changed, 67 insertions(+), 18 deletions(-)
---

diff --git a/slider.c b/slider.c

t@@ -62,7 +62,8 @@ void update_kinematics(slider s, Float dt, long int iteration)
 }
 
 // Find the relative displacement and velocity between two sliders
-void slider_displacement(slider s1, const slider s2, const Float dt, const int i)
+void slider_displacement(slider s1, const slider s2,
+        const int idx_neighbor, const int iteration)
 {
 
     // vector pointing from neighbor (s2) position to this slider position (s1)
t@@ -78,52 +79,99 @@ void slider_displacement(slider s1, const slider s2, const Float dt, const int i
     const Float3 vel_linear = subtract_float3(s1.vel, s2.vel);
 
     // relative contact interface velocity with rolling
-    const Float3 vel = add_float3(vel_linear,
+    /*const Float3 vel = add_float3(vel_linear,
             add_float3(
                 multiply_scalar_float3(dist_norm/2.,
                     cross_float3(dist, s1.angvel)),
                 multiply_scalar_float3(dist_norm/2.,
-                    cross_float3(dist, s2.angvel))));
+                    cross_float3(dist, s2.angvel))));*/
 
     // Normal component of the relative contact interface velocity
-    const Float vel_n = multiply_scalar_float3(-1.0, dot_float3(vel_linear, n));
+    const Float vel_n = -1.0*dot_float3(vel_linear, n);
 
     // Tangential component of the relative contact interface velocity
     // Hinrichsen and Wolf 2004, eq. 13.9
-    const Float3 = subtract_float3(vel, multiply_float3(n, dot_float3(n, vel)));
+    //const Float3 vel_t =
+        //subtract_float3(vel, multiply_float3_scalar(n, dot_float3(n, vel)));
 
     // read previous inter-slider distance vector
-    const Float3 dist0 = s1.neighbor_distance[i];
+    Float3 dist0;
+    if (iteration == 0)
+        dist0 = make_float3(0., 0., 0.);
+    else
+        dist0 = s1.neighbor_distance[idx_neighbor];
 
     // increment in inter-slider distance
     const Float3 ddist = subtract_float3(dist, dist0);
 
     // save current inter-slider distance
-    s1.neighbor_distance[i] = dist;
+    s1.neighbor_distance[idx_neighbor] = dist;
 
     // total relative displacement in inter-slider distance
-    s1.neighbor_relative_distance_displacement[i] =
-        add_float3(s1.neighbor_relative_distance_displacement[i], ddist);
+    if (iteration == 0)
+        s1.neighbor_relative_distance_displacement[idx_neighbor] = ddist;
+    else
+        s1.neighbor_relative_distance_displacement[idx_neighbor] =
+            add_float3(s1.neighbor_relative_distance_displacement[idx_neighbor],
+                    ddist);
 
     // total relative displacement in inter-slider distance
-    s1.neighbor_relative_distance_velocity[i] = vel_n;
-
+    s1.neighbor_relative_distance_velocity[idx_neighbor] = 
+        multiply_float3_scalar(n, vel_n);
 }
 
 
 // Resolve bond mechanics between a slider and one of its neighbors based on the 
 // incremental linear-elastic model by Potyondy and Cundall 2004
+void slider_interaction(slider s1, const slider s2, const int idx_neighbor)
+{
+
+    // determine the bond stiffness from the harmonic mean.
+    // differs from Potyondy & Stack 2004, eq. 6.
+    const Float spring_stiffness =
+        2.*s1.spring_stiffness*s2.spring_stiffness/
+        (s1.spring_stiffness + s2.spring_stiffness);
+
+
+    // bond-parallel elasticity
+    const Float3 f_n_elastic =
+        multiply_scalar_float3(bond_normal_stiffness,
+                s1.neighbor_relative_distance_displacement[idx_neighbor]);
+
+    // bond-parallel viscosity
+    const Float3 f_n_viscous =
+        multiply_scalar_float3(bond_normal_viscosity,
+                s1.neighbor_relative_distance_displacement[idx_neighbor]);
+
+    // bond-parallel force, counteracts tension and compression
+    const Float3 f_n = multiply_float3(f_n_elastic, f_n_viscous);
+
+    // add bond-parallel force to sum of forces on slider
+    s1.force = add_float3(s1.force, f_n);
+}
 
 
 // Resolve interaction between a slider and all of its neighbors
 void slider_neighbor_interaction(
         slider s,
         const slider* sliders,
-        const int N)
+        const int N,
+        const int iteration)
 {
-    int i;
-    for (i=0; i

diff --git a/slider.h b/slider.h

t@@ -30,7 +30,8 @@ typedef struct {
     Float moment_of_inertia;
 
     // inter-slider contact model parameters
-    Float spring_stiffness;  // elastic compressibility [N/m]
+    Float bond_normal_stiffness;  // Hookean elastic stiffness [N/m]
+    Float bond_normal_viscosity;  // viscosity [N/(m*s)]
 
     // slider-surface contact model parameters
     Float coulomb_friction;  // Coulomb frictional value [-]
t@@ -43,7 +44,7 @@ typedef struct {
     // relative spatial movement between this slider and its neighbors
     Float3 neighbor_distance[MAX_NEIGHBORS];
     Float3 neighbor_relative_distance_displacement[MAX_NEIGHBORS];
-    Float neighbor_relative_distance_velocity[MAX_NEIGHBORS];
+    Float3 neighbor_relative_distance_velocity[MAX_NEIGHBORS];
 
 } slider;
 

diff --git a/vector_math.h b/vector_math.h

t@@ -13,7 +13,7 @@ Float3 subtract_float3(Float3 v1, Float3 v2);
 Float3 multiply_float3(Float3 v1, Float3 v2);
 Float3 divide_float3(Float3 v1, Float3 v2);
 Float3 cross_float3(Float3 v1, Float3 v2);
-Float dot_float3(Float3 v1, Float3 v2)
+Float dot_float3(Float3 v1, Float3 v2);
 
 // vector-scalar operations
 Float3 multiply_float3_scalar(Float3 v, Float s);