tdetermine bond-parallel and bond-normal deformation inseparate functions - sli… | |
git clone git://src.adamsgaard.dk/slidergrid | |
Log | |
Files | |
Refs | |
README | |
LICENSE | |
--- | |
commit 1414f822da9c0bc70985abc65dfbca045da8cae9 | |
parent 9f1ef4954a24c00f095bfafec2afb26292b947b8 | |
Author: Anders Damsgaard <[email protected]> | |
Date: Wed, 23 Mar 2016 15:22:10 -0700 | |
determine bond-parallel and bond-normal deformation inseparate functions | |
Diffstat: | |
M Makefile | 2 +- | |
M slidergrid/slider.c | 157 +++++++++++++++++++++++++----… | |
M slidergrid/slider.h | 6 +++++- | |
3 files changed, 135 insertions(+), 30 deletions(-) | |
--- | |
diff --git a/Makefile b/Makefile | |
t@@ -10,7 +10,7 @@ ESSENTIALOBJS=$(SRCFOLDER)/main.o \ | |
$(SRCFOLDER)/simulation.o | |
BIN=example | |
-default: example | |
+default: example tests | |
example: example.o $(ESSENTIALOBJS) | |
$(CC) $(LDLIBS) $^ -o $@ | |
diff --git a/slidergrid/slider.c b/slidergrid/slider.c | |
t@@ -24,9 +24,13 @@ void initialize_slider_values(slider* s) | |
s->mass = 0.0; | |
s->moment_of_inertia = 0.0; | |
+ | |
s->bond_parallel_kv_stiffness = 0.0; | |
s->bond_parallel_kv_viscosity = 0.0; | |
+ s->bond_shear_kv_stiffness = 0.0; | |
+ s->bond_shear_kv_viscosity = 0.0; | |
+ | |
// define all entries in neighbor list as empty | |
int i; | |
for (i=0; i<MAX_NEIGHBORS; i++) { | |
t@@ -34,6 +38,7 @@ void initialize_slider_values(slider* s) | |
s->neighbor_distance[i] = zeroes_float3(); | |
s->neighbor_relative_distance_displacement[i] = zeroes_float3(); | |
s->neighbor_relative_distance_velocity[i] = zeroes_float3(); | |
+ s->neighbor_tangential_displacment[i] = zeroes_float3(); | |
} | |
} | |
t@@ -117,11 +122,9 @@ void update_kinematics(slider* s, Float dt, long int iter… | |
s->angvel = add_float3(s->angvel, dangvel_dt); | |
} | |
-// Find the relative displacement and velocity between two sliders | |
-void slider_displacement(slider* s1, const slider s2, | |
+void bond_parallel_deformation(slider* s1, const slider s2, | |
const int idx_neighbor, const int iteration) | |
{ | |
- | |
// vector pointing from neighbor (s2) position to this slider position (s1) | |
const Float3 dist = subtract_float3(s1->pos, s2.pos); | |
t@@ -135,21 +138,16 @@ void slider_displacement(slider* s1, const slider s2, | |
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 = -1.0*dot_float3(vel_linear, n); | |
- // Tangential component of the relative contact interface velocity | |
- // Hinrichsen and Wolf 2004, eq. 13.9 | |
- //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[idx_neighbor]; | |
t@@ -161,13 +159,6 @@ void slider_displacement(slider* s1, const slider s2, | |
const Float3 ddist = divide_float3_scalar( | |
subtract_float3(dist_future, dist0), 2.0); | |
- // Get displacement change from previous and current inter-slider distance | |
- // this can cause instabilities when only elasticity is included | |
- //const Float3 ddist = subtract_float3(dist, dist0); | |
- | |
- //if (iteration == 0) | |
- //ddist = zeroes_float3(); | |
- | |
// save current inter-slider distance | |
s1->neighbor_distance[idx_neighbor] = dist; | |
t@@ -183,48 +174,102 @@ void slider_displacement(slider* s1, const slider s2, | |
// total relative displacement in inter-slider distance | |
s1->neighbor_relative_distance_velocity[idx_neighbor] = | |
multiply_float3_scalar(n, vel_n); | |
+} | |
+ | |
+void bond_normal_deformation(slider* s1, const slider s2, | |
+ const int idx_neighbor, const int iteration) | |
+{ | |
+ // vector pointing from neighbor (s2) position to this slider position (s1) | |
+ const Float3 dist = subtract_float3(s1->pos, s2.pos); | |
+ | |
+ // length of inter-slider vector | |
+ const Float dist_norm = norm_float3(dist); | |
+ | |
+ // unit vector pointing from neighbor (s2) to slider s1 | |
+ const Float3 n = divide_float3_scalar(dist, dist_norm); | |
+ | |
+ // 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)))); | |
+ | |
+ // Tangential component of the relative contact interface velocity | |
+ // Hinrichsen and Wolf 2004, eq. 13.9 | |
+ const Float3 vel_t = | |
+ subtract_float3(vel, multiply_float3_scalar(n, dot_float3(n, vel))); | |
+ | |
+ // Read previous (uncorrected) tangential displacement vector | |
+ Float3 tangential_displacement0_uncor; | |
+ if (iteration == 0) | |
+ tangential_displacement0_uncor = zeroes_float3(); | |
+ else | |
+ tangential_displacement0_uncor = | |
+ s1->neighbor_tangential_displacment[idx_neighbor]; | |
+ | |
+ // Correct previous tangential displacement vector if the contact plane is | |
+ // reoriented | |
+ const Float3 tangential_displacement0 = | |
+ subtract_float3(tangential_displacement0_uncor, | |
+ multiply_float3(n, | |
+ multiply_float3(n, tangential_displacement0_uncor))); | |
- //printf("displacement = %f %f %f\n", ddist.x, ddist.y, ddist.z); | |
} | |
-// Resolve bond mechanics between a slider and one of its neighbors based on t… | |
-// incremental linear-elastic model by Potyondy and Cundall 2004 | |
-void slider_interaction(slider* s1, const slider s2, const int idx_neighbor) | |
+// Find the bond deformation | |
+void bond_deformation(slider* s1, const slider s2, | |
+ const int idx_neighbor, const int iteration) | |
{ | |
+ bond_parallel_deformation(s1, s2, idx_neighbor, iteration); | |
+ | |
+ | |
+} | |
- // determine the bond stiffness from the harmonic mean. | |
+ | |
+ | |
+ | |
+ | |
+void bond_parallel_kelvin_voigt(slider* s1, const slider s2, | |
+ const int idx_neighbor) | |
+{ | |
+ // determine the bond-parallel KV stiffness from the harmonic mean. | |
// differs from Potyondy & Stack 2004, eq. 6. | |
const Float bond_parallel_kv_stiffness = | |
2.*s1->bond_parallel_kv_stiffness*s2.bond_parallel_kv_stiffness/ | |
(s1->bond_parallel_kv_stiffness + s2.bond_parallel_kv_stiffness | |
+ 1.0e-40); | |
- // determine the bond viscosity from the harmonic mean. | |
+ // determine the bond-parallel KV viscosity from the harmonic mean. | |
const Float bond_parallel_kv_viscosity = | |
2.*s1->bond_parallel_kv_viscosity*s2.bond_parallel_kv_viscosity/ | |
(s1->bond_parallel_kv_viscosity + s2.bond_parallel_kv_viscosity | |
+ 1.0e-40); | |
- // bond-parallel elasticity | |
+ // bond-parallel Kelvin-Voigt elasticity | |
const Float3 f_n_elastic = | |
multiply_scalar_float3(bond_parallel_kv_stiffness, | |
s1->neighbor_relative_distance_displacement[idx_neighbor]); | |
- // bond-parallel viscosity | |
+ // bond-parallel Kelvin-Voigt viscosity | |
const Float3 f_n_viscous = | |
multiply_scalar_float3(bond_parallel_kv_viscosity, | |
s1->neighbor_relative_distance_velocity[idx_neighbor]); | |
- // bond-parallel force, counteracts tension and compression | |
+ // bond-parallel Kelvin-Voigt force, counteracts tension and compression | |
const Float3 f_n = multiply_scalar_float3(-1.0, | |
add_float3(f_n_elastic, f_n_viscous)); | |
- // add bond-parallel force to sum of forces on slider | |
+ // add bond-parallel Kelvin-Voigt force to sum of forces on slider | |
s1->force = add_float3(s1->force, f_n); | |
#ifdef DEBUG_SLIDER_FORCE_COMPONENTS | |
- printf("slider_interaction: f_n = %f %f %f, " | |
+ printf("slider_interaction KV-parallel: f_n = %f %f %f, " | |
"f_n_elastic = %f %f %f, f_n_viscous = %f %f %f\n", | |
f_n.x, | |
f_n.y, | |
t@@ -238,6 +283,62 @@ void slider_interaction(slider* s1, const slider s2, cons… | |
#endif | |
} | |
+void bond_shear_kelvin_voigt(slider* s1, const slider s2, | |
+ const int idx_neighbor) | |
+{ | |
+ // determine the bond-shear KV stiffness from the harmonic mean. | |
+ const Float bond_shear_kv_stiffness = | |
+ 2.*s1->bond_shear_kv_stiffness*s2.bond_shear_kv_stiffness/ | |
+ (s1->bond_shear_kv_stiffness + s2.bond_shear_kv_stiffness | |
+ + 1.0e-40); | |
+ | |
+ // determine the bond-shear KV viscosity from the harmonic mean. | |
+ const Float bond_shear_kv_viscosity = | |
+ 2.*s1->bond_shear_kv_viscosity*s2.bond_shear_kv_viscosity/ | |
+ (s1->bond_shear_kv_viscosity + s2.bond_shear_kv_viscosity | |
+ + 1.0e-40); | |
+ | |
+ // bond-parallel Kelvin-Voigt elasticity | |
+ const Float3 f_n_elastic = | |
+ multiply_scalar_float3(bond_shear_kv_stiffness, | |
+ s1->neighbor_relative_distance_displacement[idx_neighbor]); | |
+ | |
+ // bond-parallel Kelvin-Voigt viscosity | |
+ const Float3 f_n_viscous = | |
+ multiply_scalar_float3(bond_shear_kv_viscosity, | |
+ s1->neighbor_relative_distance_velocity[idx_neighbor]); | |
+ | |
+ // bond-parallel Kelvin-Voigt force, counteracts tension and compression | |
+ const Float3 f_n = multiply_scalar_float3(-1.0, | |
+ add_float3(f_n_elastic, f_n_viscous)); | |
+ | |
+ // add bond-parallel Kelvin-Voigt force to sum of forces on slider | |
+ s1->force = add_float3(s1->force, f_n); | |
+ | |
+#ifdef DEBUG_SLIDER_FORCE_COMPONENTS | |
+ printf("slider_interaction KV-shear: f_n = %f %f %f, " | |
+ "f_n_elastic = %f %f %f, f_n_viscous = %f %f %f\n", | |
+ f_n.x, | |
+ f_n.y, | |
+ f_n.z, | |
+ f_n_elastic.x, | |
+ f_n_elastic.y, | |
+ f_n_elastic.z, | |
+ f_n_viscous.x, | |
+ f_n_viscous.y, | |
+ f_n_viscous.z); | |
+#endif | |
+} | |
+ | |
+// Resolve bond mechanics between a slider and one of its neighbors based on t… | |
+// incremental linear-elastic model by Potyondy and Cundall 2004 | |
+void slider_interaction(slider* s1, const slider s2, const int idx_neighbor) | |
+{ | |
+ bond_parallel_kelvin_voigt(s1, s2, idx_neighbor); | |
+ | |
+ | |
+} | |
+ | |
// Resolve interaction between a slider and all of its neighbors | |
void slider_neighbor_interaction( | |
diff --git a/slidergrid/slider.h b/slidergrid/slider.h | |
t@@ -35,6 +35,10 @@ typedef struct { | |
Float bond_parallel_kv_stiffness; // Hookean elastic stiffness [N/m] | |
Float bond_parallel_kv_viscosity; // viscosity [N/(m*s)] | |
+ // inter-slider bond-shear Kelvin-Voigt contact model parameters | |
+ Float bond_shear_kv_stiffness; // Hookean elastic stiffness [N/m] | |
+ Float bond_shear_kv_viscosity; // viscosity [N/(m*s)] | |
+ | |
// The uniquely identifying indexes of the slider neighbors which are bond… | |
// to this slider. A value of -1 denotes an empty field. Preallocated fo… | |
// speed. | |
t@@ -44,7 +48,7 @@ typedef struct { | |
Float3 neighbor_distance[MAX_NEIGHBORS]; | |
Float3 neighbor_relative_distance_displacement[MAX_NEIGHBORS]; | |
Float3 neighbor_relative_distance_velocity[MAX_NEIGHBORS]; | |
- | |
+ Float3 neighbor_tangential_displacment[MAX_NEIGHBORS]; | |
} slider; | |