 |
tnormal-viscous contact model verified - sphere - GPU-based 3D discrete element… |
 |
git clone git://src.adamsgaard.dk/sphere |
 |
Log |
|
Files |
|
Refs |
|
LICENSE |
|
--- |
|
commit d951e69722dbdda4e924831e6ef544c468ad533e |
|
parent c62854c9bfc6be53c2155553666379bbdada700c |
|
Author: Anders Damsgaard <[email protected]> |
|
Date: Fri, 20 Jun 2014 12:34:28 +0200 |
|
|
|
normal-viscous contact model verified |
|
|
|
Diffstat: |
|
M src/contactmodels.cuh | 14 +++----------- |
|
M tests/contactmodel.py | 64 +++++++++++++++++++++++++++++… |
|
|
|
2 files changed, 67 insertions(+), 11 deletions(-) |
|
--- |
|
diff --git a/src/contactmodels.cuh b/src/contactmodels.cuh |
|
t@@ -309,9 +309,6 @@ __device__ void contactLinear(Float3* F, Float3* T, |
|
Float4 vel_b = dev_vel[idx_b_orig]; |
|
Float4 angvel4_b = dev_angvel[idx_b_orig]; |
|
|
|
- //printf("\n[%d,%d] vel = [+%e +%e]", |
|
- //idx_a_orig, idx_b_orig, vel_a, vel_b); |
|
- |
|
// Fetch previous sum of shear displacement for the contact pair |
|
Float4 delta_t0_4 = dev_delta_t[mempos]; |
|
|
|
t@@ -378,22 +375,17 @@ __device__ void contactLinear(Float3* F, Float3* T, |
|
|
|
// Normal force component: Elastic - viscous damping |
|
//f_n = (-devC_params.k_n * delta - devC_params.gamma_n * vel_n) * n; |
|
- f_n = fmax(-devC_params.k_n * delta - devC_params.gamma_n * vel_n, 0.0) * … |
|
+ f_n = fmax(0.0, -devC_params.k_n*delta + devC_params.gamma_n * vel_n) * n; |
|
+ Float f_n_length = length(f_n); |
|
|
|
// Store energy dissipated in normal viscous component |
|
// watt = gamma_n * vel_n * dx_n / dt |
|
// watt = gamma_n * vel_n * vel_n * dt / dt |
|
// watt = gamma_n * vel_n * vel_n |
|
// watt = N*m/s = N*s/m * m/s * m/s * s / s |
|
- *ev_dot += devC_params.gamma_n * vel_n * vel_n; |
|
+ *ev_dot += 0.5 * devC_params.gamma_n * vel_n * vel_n; |
|
|
|
|
|
- // Make sure the viscous damping doesn't exceed the elastic component, |
|
- // i.e. the damping factor doesn't exceed the critical damping, 2*sqrt(m*k… |
|
- if (dot(f_n, n) < 0.0f) |
|
- f_n = MAKE_FLOAT3(0.0f, 0.0f, 0.0f); |
|
- |
|
- Float f_n_length = length(f_n); |
|
|
|
// Add max. capillary force |
|
f_c = -devC_params.kappa * sqrtf(radius_a * radius_b) * n; |
|
diff --git a/tests/contactmodel.py b/tests/contactmodel.py |
|
t@@ -67,6 +67,70 @@ pytestutils.compareFloats(orig.totalKineticEnergy(), after.… |
|
"Elastic normal collision (4/4):") |
|
|
|
|
|
+## Linear viscous-elastic collisions |
|
+ |
|
+# Normal impact: Check for conservation of momentum (sum(v_i*m_i)) |
|
+orig = sphere.sim(np=2, sid='contactmodeltest') |
|
+after = sphere.sim(np=2, sid='contactmodeltest') |
|
+sphere.cleanup(orig) |
|
+orig.radius[:] = [1.0, 1.0] |
|
+orig.x[0,:] = [5.0, 5.0, 2.0] |
|
+orig.x[1,:] = [5.0, 5.0, 4.05] |
|
+v_orig = 1 |
|
+orig.vel[0,2] = v_orig |
|
+orig.defineWorldBoundaries(L=[10,10,10]) |
|
+orig.initTemporal(total = 0.1, file_dt = 0.01) |
|
+orig.gamma_n[0] = 1.0e6 |
|
+ |
|
+orig.run(verbose=False) |
|
+after.readlast(verbose=False) |
|
+#print(orig.totalKineticEnergy()) |
|
+#print(after.totalKineticEnergy()) |
|
+#print(after.totalViscousNormalEnergy()) |
|
+pytestutils.test(orig.vel[0,2] > after.vel[1,2],\ |
|
+ "Viscoelastic normal collision (1/4):") |
|
+pytestutils.compareFloats(orig.totalKineticEnergy(), |
|
+ after.totalKineticEnergy() |
|
+ + after.totalViscousNormalEnergy(), |
|
+ "Viscoelastic normal collision (2/4):", tolerance=0.05) |
|
+ |
|
+# Normal impact with different sizes: Check for conservation of momentum |
|
+orig = sphere.sim(np=2, sid='contactmodeltest') |
|
+after = sphere.sim(np=2, sid='contactmodeltest') |
|
+sphere.cleanup(orig) |
|
+orig.radius[:] = [2.0, 1.0] |
|
+orig.x[0,:] = [5.0, 5.0, 2.0] |
|
+orig.x[1,:] = [5.0, 5.0, 5.05] |
|
+orig.vel[0,2] = 1.0 |
|
+orig.defineWorldBoundaries(L=[10,10,10]) |
|
+orig.initTemporal(total = 0.1, file_dt = 0.01) |
|
+orig.gamma_n[0] = 1.0e6 |
|
+ |
|
+orig.run(verbose=False) |
|
+after.readlast(verbose=False) |
|
+pytestutils.compareFloats(orig.totalKineticEnergy(), |
|
+ after.totalKineticEnergy() |
|
+ + after.totalViscousNormalEnergy(), |
|
+ "Viscoelastic normal collision (3/4):", tolerance=0.05) |
|
+ |
|
+# Normal impact with different sizes: Check for conservation of momentum |
|
+orig = sphere.sim(np=2, sid='contactmodeltest') |
|
+after = sphere.sim(np=2, sid='contactmodeltest') |
|
+sphere.cleanup(orig) |
|
+orig.radius[:] = [1.0, 2.0] |
|
+orig.x[0,:] = [5.0, 5.0, 2.0] |
|
+orig.x[1,:] = [5.0, 5.0, 5.05] |
|
+orig.vel[0,2] = 1.0 |
|
+orig.defineWorldBoundaries(L=[10,10,10]) |
|
+orig.initTemporal(total = 0.1, file_dt = 0.01) |
|
+orig.gamma_n[0] = 1.0e6 |
|
+ |
|
+orig.run(verbose=False) |
|
+after.readlast(verbose=False) |
|
+pytestutils.compareFloats(orig.totalKineticEnergy(), |
|
+ after.totalKineticEnergy() |
|
+ + after.totalViscousNormalEnergy(), |
|
+ "Viscoelastic normal collision (4/4):", tolerance=0.05) |
|
|
|
|
|
#orig.cleanup() |
