From 066123007cd9277d4c03669707cf34ab6c05bc06 Mon Sep 17 00:00:00 2001
From: Axel Kohlmeyer <akohlmey@gmail.com>
Date: Tue, 30 May 2017 21:54:16 -0400
Subject: [PATCH] avoid undesired negative forces for high particle velocities
in granular models
---
src/GRANULAR/pair_gran_hertz_history.cpp | 2 ++
src/GRANULAR/pair_gran_hooke.cpp | 2 ++
src/GRANULAR/pair_gran_hooke_history.cpp | 2 ++
3 files changed, 6 insertions(+)
diff --git a/src/GRANULAR/pair_gran_hertz_history.cpp b/src/GRANULAR/pair_gran_hertz_history.cpp
index e52aac10db..14cbc8eb0a 100644
--- a/src/GRANULAR/pair_gran_hertz_history.cpp
+++ b/src/GRANULAR/pair_gran_hertz_history.cpp
@@ -183,6 +183,7 @@ void PairGranHertzHistory::compute(int eflag, int vflag)
ccel = kn*(radsum-r)*rinv - damp;
polyhertz = sqrt((radsum-r)*radi*radj / radsum);
ccel *= polyhertz;
+ if (ccel < 0.0) ccel = 0.0;
// relative velocities
@@ -390,6 +391,7 @@ double PairGranHertzHistory::single(int i, int j, int itype, int jtype,
ccel = kn*(radsum-r)*rinv - damp;
polyhertz = sqrt((radsum-r)*radi*radj / radsum);
ccel *= polyhertz;
+ if (ccel < 0.0) ccel = 0.0;
// relative velocities
diff --git a/src/GRANULAR/pair_gran_hooke.cpp b/src/GRANULAR/pair_gran_hooke.cpp
index 66e1c3dd7a..2d6d113ab7 100644
--- a/src/GRANULAR/pair_gran_hooke.cpp
+++ b/src/GRANULAR/pair_gran_hooke.cpp
@@ -161,6 +161,7 @@ void PairGranHooke::compute(int eflag, int vflag)
damp = meff*gamman*vnnr*rsqinv;
ccel = kn*(radsum-r)*rinv - damp;
+ if (ccel < 0.0) ccel = 0.0;
// relative velocities
@@ -302,6 +303,7 @@ double PairGranHooke::single(int i, int j, int itype, int jtype, double rsq,
damp = meff*gamman*vnnr*rsqinv;
ccel = kn*(radsum-r)*rinv - damp;
+ if (ccel < 0.0) ccel = 0.0;
// relative velocities
diff --git a/src/GRANULAR/pair_gran_hooke_history.cpp b/src/GRANULAR/pair_gran_hooke_history.cpp
index e9662c9e73..1378709cb5 100644
--- a/src/GRANULAR/pair_gran_hooke_history.cpp
+++ b/src/GRANULAR/pair_gran_hooke_history.cpp
@@ -223,6 +223,7 @@ void PairGranHookeHistory::compute(int eflag, int vflag)
damp = meff*gamman*vnnr*rsqinv;
ccel = kn*(radsum-r)*rinv - damp;
+ if (ccel < 0.0) ccel = 0.0;
// relative velocities
@@ -687,6 +688,7 @@ double PairGranHookeHistory::single(int i, int j, int itype, int jtype,
damp = meff*gamman*vnnr*rsqinv;
ccel = kn*(radsum-r)*rinv - damp;
+ if (ccel < 0.0) ccel = 0.0;
// relative velocities
--
GitLab