diff options
Diffstat (limited to 'src/particles/forces/hertzMindlinDeresiewiczForCombWithCollisionModel3D.hh')
| -rwxr-xr-x | src/particles/forces/hertzMindlinDeresiewiczForCombWithCollisionModel3D.hh | 235 |
1 files changed, 235 insertions, 0 deletions
diff --git a/src/particles/forces/hertzMindlinDeresiewiczForCombWithCollisionModel3D.hh b/src/particles/forces/hertzMindlinDeresiewiczForCombWithCollisionModel3D.hh new file mode 100755 index 0000000..eba2f31 --- /dev/null +++ b/src/particles/forces/hertzMindlinDeresiewiczForCombWithCollisionModel3D.hh @@ -0,0 +1,235 @@ +/* + * 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, template< + typename W> class 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, template< + typename W> class 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, template< + typename W> class DESCRIPTOR> +void HertzMindlinDeresiewicz3D<T, PARTICLETYPE, DESCRIPTOR>::computeForce( + typename std::deque<PARTICLETYPE<T> >::iterator p, int pInt, + ParticleSystem3D<T, PARTICLETYPE>& pSys, T force[3]) +{ + if (p->getSActivity() > 1) { + + 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); + + /// Limit Overlap + // T deltaMax = 0.03 * (p2->getRad() + p->getRad()) ; + // if (delta > deltaMax ) { + + // T dpos[3] = {T(0), T(0), T(0) } ; + + // for (int i = 0; i <= 2; i++) { + + // dpos[i] = _normal[i] * 0.5 * (delta - deltaMax) ; + // p->getPos()[i] -= 1.* dpos[i]; + // p2->getPos()[i] += 1.* dpos[i]; + // } + // delta = deltaMax * (p2->getRad() + p->getRad()); + // } + + // 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]); // gehört das Minus hier hin? + _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 ; + + if ((p->getSActivity() || p2->getSActivity()) == 3) { + p->setSActivity(3); + p2->setSActivity(3); + } + if ((p->getSActivity() == 4) && (p2->getSActivity() != (4 || 3))) { + p2->setSActivity(3); + } + if ((p2->getSActivity() == 4) && (p->getSActivity() != (4 || 3))) { + p->setSActivity(3); + } + } + } + } + } +} + +#endif + + |