Initialize at openlb-1-3

1 files changed, 202 insertions, 0 deletions

diff --git a/src/particles/forces/hertzMindlinDeresiewicz3D.hh b/src/particles/forces/hertzMindlinDeresiewicz3D.hh
new file mode 100755
index 0000000..78bcf2e
--- /dev/null
+++ b/src/particles/forces/hertzMindlinDeresiewicz3D.hh
@@ -0,0 +1,202 @@
+/*
+ *  Copyright (C) 2015 Marie-Luise Maier, Mathias J. Krause, Sascha Janz
+ *  E-mail contact: info@openlb.net
+ *  The most recent release of OpenLB can be downloaded at
+ *  <http://www.openlb.net/>
+ *
+ *  This program is free software; you can redistribute it and/or
+ *  modify it under the terms of the GNU General Public License
+ *  as published by the Free Software Foundation; either version 2
+ *  of the License, or (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public
+ *  License along with this program; if not, write to the Free
+ *  Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+ *  Boston, MA  02110-1301, USA.
+ */
+
+/** Alberto Di Renzo, Francesco Paolo Di Maio:
+ * "Comparison of contact-force models for the simulation of collisions in
+ * DEM-based granular ow codes",
+ * Chemical Engineering Science 59 (2004) 525 - 541
+ */
+
+#ifndef HERTZMINDLINDERESIEWICZ3D_HH
+#define HERTZMINDLINDERESIEWICZ3D_HH
+
+#include <cmath>
+#include "hertzMindlinDeresiewicz3D.h"
+
+namespace olb {
+
+template<typename T, template<typename U> class PARTICLETYPE, typename DESCRIPTOR>
+HertzMindlinDeresiewicz3D<T, PARTICLETYPE, DESCRIPTOR>::HertzMindlinDeresiewicz3D(
+  T G1, T G2, T v1, T v2, T scale1, T scale2, bool validationKruggelEmden) :
+  Force3D<T, PARTICLETYPE>(), _G1(G1), _G2(G2), _v1(v1), _v2(v2), _scale1(
+    scale1), _scale2(scale2), _validationKruggelEmden(validationKruggelEmden)
+{
+  // E-Modul Particle
+  E1 = 2 * (1 + _v1) * _G1;
+  E2 = 2 * (1 + _v2) * _G2;
+
+  // equivalent combined E-Modul
+  eE = (1 - pow(_v1, 2)) / E1 + (1 - pow(_v2, 2)) / E2;
+  eE = 1 / eE;
+
+  // equivalent combined E-Modul
+  eG = (2.0 - _v1) / _G1 + (2 - _v2) / _G2;
+  eG = 1. / eG;
+}
+
+template<typename T, template<typename U> class PARTICLETYPE, typename DESCRIPTOR>
+void HertzMindlinDeresiewicz3D<T, PARTICLETYPE, DESCRIPTOR>::applyForce(
+  typename std::deque<PARTICLETYPE<T> >::iterator p, int pInt,
+  ParticleSystem3D<T, PARTICLETYPE>& pSys)
+{
+  T force[3] = {T(), T(), T()};
+
+  computeForce(p, pInt, pSys, force);
+}
+
+template<typename T, template<typename U> class PARTICLETYPE, typename DESCRIPTOR>
+void HertzMindlinDeresiewicz3D<T, PARTICLETYPE, DESCRIPTOR>::computeForce(
+  typename std::deque<PARTICLETYPE<T> >::iterator p, int pInt,
+  ParticleSystem3D<T, PARTICLETYPE>& pSys, T force[3])
+{
+
+  std::vector<std::pair<size_t, T>> ret_matches;
+  // kind of contactDetection has to be chosen in application
+  pSys.getContactDetection()->getMatches(pInt, ret_matches);
+
+  PARTICLETYPE<T>* p2 = NULL;
+
+  // iterator walks through number of neighbored particles = ret_matches
+  for (const auto& it : ret_matches) {
+
+    if (!util::nearZero(it.second)) {
+
+      p2 = &pSys[it.first];
+
+      // overlap
+      T delta = (p2->getRad() + p->getRad()) - sqrt(it.second);
+
+      // equivalent mass
+      T M = p->getMass() * p2->getMass() / (p->getMass() + p2->getMass());
+      // equivalent radius
+      T R = p->getRad() * p2->getRad() / (p->getRad() + p2->getRad());
+      // relative velocity
+      std::vector < T > _velR(3, T());
+      _velR[0] = -(p2->getVel()[0] - p->getVel()[0]);
+      _velR[1] = -(p2->getVel()[1] - p->getVel()[1]);
+      _velR[2] = -(p2->getVel()[2] - p->getVel()[2]);
+
+      std::vector < T > _d(3, T());
+      std::vector < T > _normal(3, T());
+
+      //_d: vector from particle1 to particle2
+      _d[0] = p2->getPos()[0] - p->getPos()[0];
+      _d[1] = p2->getPos()[1] - p->getPos()[1];
+      _d[2] = p2->getPos()[2] - p->getPos()[2];
+
+      if ( !util::nearZero(util::norm(_d)) ) {
+        _normal = util::normalize(_d);
+      }
+      else {
+        return;
+      }
+
+      Vector<T, 3> d_(_d);
+      Vector<T, 3> velR_(_velR);
+      T dot = velR_[0] * _normal[0] + velR_[1] * _normal[1] + velR_[2] * _normal[2];
+
+      // normal part of relative velocity
+      // normal relative to surface of particles at contact point
+      std::vector < T > _velN(3, T());
+      _velN[0] = dot * _normal[0];
+      _velN[1] = dot * _normal[1];
+      _velN[2] = dot * _normal[2];
+
+      // tangential part of relative velocity
+      // tangential relative to surface of particles at contact point
+      std::vector < T > _velT(3, T());
+      _velT[0] = _velR[0] - _velN[0];
+      _velT[1] = _velR[1] - _velN[1];
+      _velT[2] = _velR[2] - _velN[2];
+
+      if (delta > 0.) {
+
+        // Force normal
+        // spring constant in normal direction
+        // (Alberto Di Renzo, Francesco Paolo Di Maio, Chemical Engineering Science 59 (2004) 525 - 541)
+        // constant kn from H. Kruggel-Endem
+        T kn = 4 / 3. * sqrt(R) * eE;
+        if (_validationKruggelEmden) {
+          kn = 7.35e9;  // to compare to Kruggel-Emden
+        }
+
+        // part of mechanical force of spring in normal direction
+        // Hertz Contact (P. A. Langston, Powder Technology 85 (1995))
+        std::vector < T > Fs_n(3, T());
+        Fs_n[0] = -kn * pow(delta, 1.5) * _normal[0];
+        Fs_n[1] = -kn * pow(delta, 1.5) * _normal[1];
+        Fs_n[2] = -kn * pow(delta, 1.5) * _normal[2];
+
+        // part of mechanical force of damper in normal direction
+        // damped linear spring (Cundall, Strack 1979)
+        // (K.W. Chu, A.B. Yu, Powder Technology 179 (2008) 104 – 114)
+        // damper constant in normal direction
+        // constant eta_n from H. Kruggel-Endem
+        T eta_n = 0.3 * sqrt(4.5 * M * sqrt(delta) * kn);
+        if (_validationKruggelEmden) {
+          eta_n = 1.96e5;  // to compare to Kruggel-Emden
+        }
+
+        std::vector < T > Fd_n(3, T());
+        Fd_n[0] = -eta_n * _velN[0] * sqrt(delta);
+        Fd_n[1] = -eta_n * _velN[1] * sqrt(delta);
+        Fd_n[2] = -eta_n * _velN[2] * sqrt(delta);
+
+        std::vector < T > F_n(3, T());
+        F_n[0] = Fs_n[0] + Fd_n[0];
+        F_n[1] = Fs_n[1] + Fd_n[1];
+        F_n[2] = Fs_n[2] + Fd_n[2];
+
+        // Force tangential
+        // spring constant in tangential direction
+        // (N.G. Deen, Chemical Engineering Science 62 (2007) 28 - 44)
+        T kt = 2 * sqrt(2 * R) * _G1 / (2 - _v1) * pow(delta, 0.5);
+
+        // damper constant in normal direction
+        T eta_t = 2 * sqrt(2. / 7. * M * kt);
+
+        // part of mechanical force of damper in tangential direction
+        std::vector < T > F_t(3, T());
+        F_t[0] = -eta_t * _velT[0];
+        F_t[1] = -eta_t * _velT[1];
+        F_t[2] = -eta_t * _velT[2];
+
+        // entire force
+        // factor _scale to prevent instability
+        force[0] = _scale1 * F_n[0] + _scale2 * F_t[0];
+        force[1] = _scale1 * F_n[1] + _scale2 * F_t[1];
+        force[2] = _scale1 * F_n[2] + _scale2 * F_t[2];
+
+        p->getForce()[0] += force[0] * 0.5 ;
+        p->getForce()[1] += force[1] * 0.5 ;
+        p->getForce()[2] += force[2] * 0.5 ;
+        p2->getForce()[0] -= force[0] * 0.5 ;
+        p2->getForce()[1] -= force[1] * 0.5 ;
+        p2->getForce()[2] -= force[2] * 0.5 ;
+      }
+    }
+  }
+}
+
+}
+
+#endif