commit c1d935d91f4cf966769161facb856ed8c426f859
parent 25a1b1a29874237888551abd9d157460aedac79f
Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Wed, 30 Mar 2016 10:51:54 -0700
add viscous and non-viscous damping
Diffstat:
4 files changed, 71 insertions(+), 17 deletions(-)
diff --git a/slidergrid/main.c b/slidergrid/main.c
@@ -175,6 +175,18 @@ int main(int argc, char** argv)
#endif
}
+ // add global viscous damping
+ for (i=0; i<sim.N; i++)
+ if (sim.sliders[i].damping_viscosity_linear > 1.0e-16 ||
+ sim.sliders[i].damping_viscosity_angular > 1.0e-16)
+ slider_global_viscous_damping(&sim.sliders[i]);
+
+ // add global non-viscous damping
+ for (i=0; i<sim.N; i++)
+ if (sim.sliders[i].damping_coefficient > 1.0e-16)
+ slider_nonviscous_damping(&sim.sliders[i]);
+
+
for (i=0; i<sim.N; i++)
update_kinematics(&sim.sliders[i], sim.dt, sim.iteration);
diff --git a/slidergrid/simulation.c b/slidergrid/simulation.c
@@ -220,7 +220,7 @@ int save_sliders_to_vtk_file(
// angular positions
fprintf(f, " <DataArray NumberOfComponents=\"3\" type=\"Float32\" "
- "Name=\"Angular position [m]\" format=\"ascii\">\n");
+ "Name=\"Angular position [rad]\" format=\"ascii\">\n");
fprintf(f, " ");
for (i=0; i<N; i++)
fprintf(f, "%f %f %f ",
diff --git a/slidergrid/slider.c b/slidergrid/slider.c
@@ -1,4 +1,5 @@
#include <stdio.h>
+#include <math.h>
#include "typedefs.h"
#include "slider.h"
#include "vector_math.h"
@@ -31,6 +32,10 @@ void initialize_slider_values(slider* s)
s->bond_shear_kv_stiffness = 0.0;
s->bond_shear_kv_viscosity = 0.0;
+ s->damping_viscosity_linear = 0.0;
+ s->damping_viscosity_angular = 0.0;
+ s->damping_coefficient = 0.0;
+
// define all entries in neighbor list as empty
int i;
for (i=0; i<MAX_NEIGHBORS; i++) {
@@ -137,14 +142,6 @@ void bond_parallel_deformation(slider* s1, const slider s2,
// relative contact interface velocity w/o rolling
const Float3 vel_linear = subtract_float3(s1->vel, s2.vel);
- // relative contact interface velocity with rolling
- 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))));
-
// Normal component of the relative contact interface velocity
const Float vel_n = -1.0*dot_float3(vel_linear, n);
@@ -219,6 +216,15 @@ void bond_normal_deformation(slider* s1, const slider s2,
multiply_float3(n,
multiply_float3(n, tangential_displacement0_uncor)));
+ // project future tangential displacement
+
+
+ // determine dtangential_displacement by central differences
+
+
+ // use dtangential_displacement for elastic response
+
+ // use vel_t for viscous response
}
@@ -227,14 +233,10 @@ void bond_deformation(slider* s1, const slider s2,
const int idx_neighbor, const int iteration)
{
bond_parallel_deformation(s1, s2, idx_neighbor, iteration);
-
-
+ //bond_normal_deformation(s1, s2, idx_neighbor, iteration);
}
-
-
-
void bond_parallel_kelvin_voigt(slider* s1, const slider s2,
const int idx_neighbor)
{
@@ -335,8 +337,7 @@ void bond_shear_kelvin_voigt(slider* s1, const slider s2,
void slider_interaction(slider* s1, const slider s2, const int idx_neighbor)
{
bond_parallel_kelvin_voigt(s1, s2, idx_neighbor);
-
-
+ //bond_normal_kelvin_voigt(s1, s2, idx_neighbor);
}
@@ -360,7 +361,7 @@ void slider_neighbor_interaction(
if (s->neighbors[idx_neighbor] != -1) {
- slider_displacement(
+ bond_deformation(
s, sliders[s->neighbors[idx_neighbor]],
idx_neighbor, iteration);
@@ -379,3 +380,36 @@ void slider_neighbor_interaction(
}
}
}
+
+// add frequency-independent viscosity to sum of forces and torques to damp
+// reflected waves from the edge of the lattice
+// Mora and Place (1994), Place et al (2002)
+void slider_viscous_damping(slider* s)
+{
+ s->force.x += - s->damping_viscosity_linear*s->vel.x;
+ s->force.y += - s->damping_viscosity_linear*s->vel.y;
+ s->force.z += - s->damping_viscosity_linear*s->vel.z;
+
+ s->torque.x += - s->damping_viscosity_angular*s->vel.x;
+ s->torque.y += - s->damping_viscosity_angular*s->vel.y;
+ s->torque.z += - s->damping_viscosity_angular*s->vel.z;
+}
+
+// add local non-viscous damping, Potyondy and Cundall (2004)
+void slider_nonviscous_damping(slider* s)
+{
+ s->force.x +=
+ - s->damping_coefficient*fabs(s->force.x)*copysign(1.0, s->vel.x);
+ s->force.y +=
+ - s->damping_coefficient*fabs(s->force.y)*copysign(1.0, s->vel.y);
+ s->force.z +=
+ - s->damping_coefficient*fabs(s->force.z)*copysign(1.0, s->vel.z);
+
+ s->torque.x +=
+ - s->damping_coefficient*fabs(s->torque.x)*copysign(1.0, s->angvel.x);
+ s->torque.y +=
+ - s->damping_coefficient*fabs(s->torque.y)*copysign(1.0, s->angvel.y);
+ s->torque.z +=
+ - s->damping_coefficient*fabs(s->torque.z)*copysign(1.0, s->angvel.z);
+
+}
diff --git a/slidergrid/slider.h b/slidergrid/slider.h
@@ -39,6 +39,11 @@ typedef struct {
Float bond_shear_kv_stiffness; // Hookean elastic stiffness [N/m]
Float bond_shear_kv_viscosity; // viscosity [N/(m*s)]
+ // Damping parameters
+ Float damping_viscosity_linear; // Linear velocity damping [N/(m/s)]
+ Float damping_viscosity_angular; // Angular velocity damping [N*m/(rad/s)]
+ Float damping_coefficient; // Dimensionless damping [-]
+
// The uniquely identifying indexes of the slider neighbors which are bonded
// to this slider. A value of -1 denotes an empty field. Preallocated for
// speed.
@@ -65,4 +70,7 @@ void slider_neighbor_interaction(
const int N,
const int iteration);
+void slider_viscous_damping(slider* s);
+void slider_nonviscous_damping(slider* s);
+
#endif
