From 94d3e79a8617f88dc0219cfdeedfa3147833719d Mon Sep 17 00:00:00 2001 From: Adrian Kummerlaender Date: Mon, 24 Jun 2019 14:43:36 +0200 Subject: Initialize at openlb-1-3 --- src/particles/forces/MakeHeader | 36 ++++ src/particles/forces/buoyancyForce3D.cpp | 35 +++ src/particles/forces/buoyancyForce3D.h | 58 +++++ src/particles/forces/buoyancyForce3D.hh | 67 ++++++ src/particles/forces/force3D.cpp | 34 +++ src/particles/forces/force3D.h | 58 +++++ src/particles/forces/force3D.hh | 50 +++++ src/particles/forces/forceFromExtField3D.cpp | 35 +++ src/particles/forces/forceFromExtField3D.h | 57 +++++ src/particles/forces/forceFromExtField3D.hh | 66 ++++++ src/particles/forces/forces.h | 42 ++++ src/particles/forces/forces.hh | 40 ++++ src/particles/forces/hertzMindlinDeresiewicz3D.cpp | 36 ++++ src/particles/forces/hertzMindlinDeresiewicz3D.h | 75 +++++++ src/particles/forces/hertzMindlinDeresiewicz3D.hh | 202 ++++++++++++++++++ ...ndlinDeresiewiczForCombWithCollisionModel3D.cpp | 36 ++++ ...MindlinDeresiewiczForCombWithCollisionModel3D.h | 76 +++++++ ...indlinDeresiewiczForCombWithCollisionModel3D.hh | 235 +++++++++++++++++++++ src/particles/forces/interpMagForceForMagP3D.cpp | 40 ++++ src/particles/forces/interpMagForceForMagP3D.h | 59 ++++++ src/particles/forces/interpMagForceForMagP3D.hh | 158 ++++++++++++++ src/particles/forces/linearContactForce3D.cpp | 36 ++++ src/particles/forces/linearContactForce3D.h | 76 +++++++ src/particles/forces/linearContactForce3D.hh | 208 ++++++++++++++++++ .../linearDampingForceForMagDipoleMoment3D.cpp | 37 ++++ .../linearDampingForceForMagDipoleMoment3D.h | 55 +++++ .../linearDampingForceForMagDipoleMoment3D.hh | 63 ++++++ src/particles/forces/magneticForceForMagP3D.cpp | 35 +++ src/particles/forces/magneticForceForMagP3D.h | 64 ++++++ src/particles/forces/magneticForceForMagP3D.hh | 85 ++++++++ src/particles/forces/magneticForceFromHField3D.cpp | 35 +++ src/particles/forces/magneticForceFromHField3D.h | 57 +++++ src/particles/forces/magneticForceFromHField3D.hh | 119 +++++++++++ src/particles/forces/module.mk | 27 +++ src/particles/forces/stokesDragForce3D.cpp | 36 ++++ src/particles/forces/stokesDragForce3D.h | 71 +++++++ src/particles/forces/stokesDragForce3D.hh | 127 +++++++++++ .../forces/stokesDragForceForHomVelField3D.cpp | 38 ++++ .../forces/stokesDragForceForHomVelField3D.h | 54 +++++ .../forces/stokesDragForceForHomVelField3D.hh | 57 +++++ src/particles/forces/transferExternalForce3D.cpp | 33 +++ src/particles/forces/transferExternalForce3D.h | 48 +++++ src/particles/forces/transferExternalForce3D.hh | 54 +++++ src/particles/forces/weightForce3D.cpp | 34 +++ src/particles/forces/weightForce3D.h | 50 +++++ src/particles/forces/weightForce3D.hh | 50 +++++ 46 files changed, 3044 insertions(+) create mode 100644 src/particles/forces/MakeHeader create mode 100644 src/particles/forces/buoyancyForce3D.cpp create mode 100644 src/particles/forces/buoyancyForce3D.h create mode 100644 src/particles/forces/buoyancyForce3D.hh create mode 100644 src/particles/forces/force3D.cpp create mode 100644 src/particles/forces/force3D.h create mode 100644 src/particles/forces/force3D.hh create mode 100644 src/particles/forces/forceFromExtField3D.cpp create mode 100755 src/particles/forces/forceFromExtField3D.h create mode 100755 src/particles/forces/forceFromExtField3D.hh create mode 100644 src/particles/forces/forces.h create mode 100644 src/particles/forces/forces.hh create mode 100644 src/particles/forces/hertzMindlinDeresiewicz3D.cpp create mode 100755 src/particles/forces/hertzMindlinDeresiewicz3D.h create mode 100755 src/particles/forces/hertzMindlinDeresiewicz3D.hh create mode 100644 src/particles/forces/hertzMindlinDeresiewiczForCombWithCollisionModel3D.cpp create mode 100755 src/particles/forces/hertzMindlinDeresiewiczForCombWithCollisionModel3D.h create mode 100755 src/particles/forces/hertzMindlinDeresiewiczForCombWithCollisionModel3D.hh create mode 100644 src/particles/forces/interpMagForceForMagP3D.cpp create mode 100644 src/particles/forces/interpMagForceForMagP3D.h create mode 100644 src/particles/forces/interpMagForceForMagP3D.hh create mode 100644 src/particles/forces/linearContactForce3D.cpp create mode 100755 src/particles/forces/linearContactForce3D.h create mode 100755 src/particles/forces/linearContactForce3D.hh create mode 100644 src/particles/forces/linearDampingForceForMagDipoleMoment3D.cpp create mode 100644 src/particles/forces/linearDampingForceForMagDipoleMoment3D.h create mode 100644 src/particles/forces/linearDampingForceForMagDipoleMoment3D.hh create mode 100644 src/particles/forces/magneticForceForMagP3D.cpp create mode 100755 src/particles/forces/magneticForceForMagP3D.h create mode 100755 src/particles/forces/magneticForceForMagP3D.hh create mode 100644 src/particles/forces/magneticForceFromHField3D.cpp create mode 100755 src/particles/forces/magneticForceFromHField3D.h create mode 100755 src/particles/forces/magneticForceFromHField3D.hh create mode 100644 src/particles/forces/module.mk create mode 100644 src/particles/forces/stokesDragForce3D.cpp create mode 100644 src/particles/forces/stokesDragForce3D.h create mode 100644 src/particles/forces/stokesDragForce3D.hh create mode 100644 src/particles/forces/stokesDragForceForHomVelField3D.cpp create mode 100644 src/particles/forces/stokesDragForceForHomVelField3D.h create mode 100644 src/particles/forces/stokesDragForceForHomVelField3D.hh create mode 100644 src/particles/forces/transferExternalForce3D.cpp create mode 100644 src/particles/forces/transferExternalForce3D.h create mode 100644 src/particles/forces/transferExternalForce3D.hh create mode 100644 src/particles/forces/weightForce3D.cpp create mode 100644 src/particles/forces/weightForce3D.h create mode 100644 src/particles/forces/weightForce3D.hh (limited to 'src/particles/forces') diff --git a/src/particles/forces/MakeHeader b/src/particles/forces/MakeHeader new file mode 100644 index 0000000..96b8e6b --- /dev/null +++ b/src/particles/forces/MakeHeader @@ -0,0 +1,36 @@ +# This file is part of the OpenLB library +# +# Copyright (C) 2007 Mathias Krause +# E-mail contact: info@openlb.net +# The most recent release of OpenLB can be downloaded at +# +# +# 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. + + +generic := + +precompiled := force3D \ + stokesDragForce3D \ + weightForce3D \ + buoyancyForce3D \ + hertzMindlinDeresiewicz3D \ + transferExternalForce3D \ + magneticForceFromHField3D \ + linearDampingForceForMagDipoleMoment3D \ + magneticForceForMagP3D \ + interpMagForceForMagP3D + diff --git a/src/particles/forces/buoyancyForce3D.cpp b/src/particles/forces/buoyancyForce3D.cpp new file mode 100644 index 0000000..450597a --- /dev/null +++ b/src/particles/forces/buoyancyForce3D.cpp @@ -0,0 +1,35 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2016 Mathias J. Krause, Marie-Luise Maier + * E-mail contact: info@openlb.net + * The most recent release of OpenLB can be downloaded at + * + * + * 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. +*/ + +#include "dynamics/latticeDescriptors.h" + +#include "particles/particle3D.h" +#include "particles/particle3D.hh" +#include "buoyancyForce3D.h" +#include "buoyancyForce3D.hh" + +namespace olb { + +template class BuoyancyForce3D>; + +} // namespace olb diff --git a/src/particles/forces/buoyancyForce3D.h b/src/particles/forces/buoyancyForce3D.h new file mode 100644 index 0000000..27bb41a --- /dev/null +++ b/src/particles/forces/buoyancyForce3D.h @@ -0,0 +1,58 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2013 Thomas Henn, Mathias J. Krause + * E-mail contact: info@openlb.net + * The most recent release of OpenLB can be downloaded at + * + * + * 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. + */ + +#ifndef BUOYANCYFORCE3D_H_ +#define BUOYANCYFORCE3D_H_ + +#include "particles/particleSystem3D.h" +#include "forces.h" +#include "core/unitConverter.h" +//#include "../../functors/superLatticeLocalF3D.h" + + +namespace olb { + +template class PARTICLETYPE> +class ParticleSystem3D; + +template class PARTICLETYPE, typename DESCRIPTOR> +class BuoyancyForce3D: public Force3D { + +public: + BuoyancyForce3D(UnitConverter const& converter, std::vector direction, + T g = 9.81); + ~BuoyancyForce3D() override { }; + void applyForce(typename std::deque >::iterator p, int pInt, + ParticleSystem3D& psSys) override; + +// void computeForce(int pInt, ParticleSystem3D* psSys, +// T force[3]); +private: + std::vector _direction; + T _g; + T _physDensity; +}; + +} + +#endif /* BUOYANCYFORCE3D_H_ */ diff --git a/src/particles/forces/buoyancyForce3D.hh b/src/particles/forces/buoyancyForce3D.hh new file mode 100644 index 0000000..b05dada --- /dev/null +++ b/src/particles/forces/buoyancyForce3D.hh @@ -0,0 +1,67 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2016 Thomas Henn + * E-mail contact: info@openlb.net + * The most recent release of OpenLB can be downloaded at + * + * + * 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. + */ + +#ifndef BUOYANCYFORCE_3D_HH +#define BUOYANCYFORCE_3D_HH + +#include +#include "buoyancyForce3D.h" + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +namespace olb { + +template class PARTICLETYPE, typename DESCRIPTOR> +BuoyancyForce3D::BuoyancyForce3D( + UnitConverter const& converter, + std::vector direction, T g) : + Force3D(), + _direction(direction), _g(g), _physDensity(converter.getPhysDensity() ) +{ + T directionNorm = sqrt( + pow(_direction[0], 2.) + pow(_direction[1], 2.) + + pow(_direction[2], 2.)); + for (int i = 0; i < 3; ++i) { + _direction[i] /= directionNorm; + } +} + +template class PARTICLETYPE, typename DESCRIPTOR> +void BuoyancyForce3D::applyForce( + typename std::deque >::iterator p, int pInt, + ParticleSystem3D& psSys) +{ + // weight force of fluid in similar volume like particle that replaces the fluid + // acts in direction opposite to weight force of particle + + T factor = 4. / 3. * M_PI * std::pow(p->getRad(), 3) * _g + * _physDensity; + for (int j = 0; j < 3; ++j) { + p->getForce()[j] -= factor * _direction[j]; + } +} + +} +#endif /* BUOYANCYFORCE_3D_HH */ diff --git a/src/particles/forces/force3D.cpp b/src/particles/forces/force3D.cpp new file mode 100644 index 0000000..e716730 --- /dev/null +++ b/src/particles/forces/force3D.cpp @@ -0,0 +1,34 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2016 Mathias J. Krause, Marie-Luise Maier + * E-mail contact: info@openlb.net + * The most recent release of OpenLB can be downloaded at + * + * + * 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. +*/ + +#include"particles/particle3D.h" +#include"particles/particle3D.hh" +#include"force3D.h" +#include"force3D.hh" + +namespace olb { + +template class Force3D; +template class Force3D; + +} // namespace olb diff --git a/src/particles/forces/force3D.h b/src/particles/forces/force3D.h new file mode 100644 index 0000000..4073363 --- /dev/null +++ b/src/particles/forces/force3D.h @@ -0,0 +1,58 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2016 Thomas Henn + * E-mail contact: info@openlb.net + * The most recent release of OpenLB can be downloaded at + * + * + * 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. +*/ + + +#ifndef FORCE_3D_H +#define FORCE_3D_H + +#include +#include "io/ostreamManager.h" +#include "particles/particleSystem3D.h" + +namespace olb { + +template class PARTICLETYPE> +class ParticleSystem3D; + + +/** + * Prototype for all particle forces + */ + +template class PARTICLETYPE> +class Force3D { +public: + Force3D(); + Force3D(Force3D&); + Force3D(const Force3D&); + + virtual ~Force3D() {}; + virtual void applyForce(typename std::deque >::iterator p, int pInt, ParticleSystem3D& psSys)=0; + +protected: + mutable OstreamManager clout; + +}; + +} +#endif /* FORCE3D_H */ diff --git a/src/particles/forces/force3D.hh b/src/particles/forces/force3D.hh new file mode 100644 index 0000000..0e27d25 --- /dev/null +++ b/src/particles/forces/force3D.hh @@ -0,0 +1,50 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2016 Thomas Henn + * E-mail contact: info@openlb.net + * The most recent release of OpenLB can be downloaded at + * + * + * 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. +*/ + + +#ifndef FORCE_3D_HH +#define FORCE_3D_HH + +#include "force3D.h" + +namespace olb { + +template class PARTICLETYPE> +Force3D::Force3D(Force3D& f) : + clout(f.clout) +{ +} + +template class PARTICLETYPE> +Force3D::Force3D(const Force3D& f) : + clout(f.clout) +{ +} + +template class PARTICLETYPE> +Force3D::Force3D(): clout(std::cout,"Force3D") +{ +} + +} +#endif /* FORCE3D_HH */ diff --git a/src/particles/forces/forceFromExtField3D.cpp b/src/particles/forces/forceFromExtField3D.cpp new file mode 100644 index 0000000..7bd7d20 --- /dev/null +++ b/src/particles/forces/forceFromExtField3D.cpp @@ -0,0 +1,35 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2018 Marie-Luise Maier + * E-mail contact: info@openlb.net + * The most recent release of OpenLB can be downloaded at + * + * + * 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. +*/ + +#include "dynamics/latticeDescriptors.h" + +#include "particles/particle3D.h" +#include "particles/particle3D.hh" +#include "forceFromExtField3D.h" +#include "forceFromExtField3D.hh" + +namespace olb { + +template class ForceFromExtField3D>; + +} // namespace olb diff --git a/src/particles/forces/forceFromExtField3D.h b/src/particles/forces/forceFromExtField3D.h new file mode 100755 index 0000000..20cbaa6 --- /dev/null +++ b/src/particles/forces/forceFromExtField3D.h @@ -0,0 +1,57 @@ +/* + * Copyright (C) 2018 Marie-Luise Maier + * E-mail contact: info@openlb.net + * The most recent release of OpenLB can be downloaded at + * + * + * 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. + */ + +#ifndef ForceFromExtField3D_H +#define ForceFromExtField3D_H + +#include "functors/lattice/superLatticeLocalF3D.h" +#include "particles/particleSystem3D.h" +#include "force3D.h" + +namespace olb { + +template class PARTICLETYPE> +class ParticleSystem3D; + +template class PARTICLETYPE, typename DESCRIPTOR> +class ForceFromExtField3D: public Force3D { + +public: + // todo rename in ForceFromFunction3D or create new force + ForceFromExtField3D( +// SuperLattice3D& sLattice, + //SuperLatticeField3D& sLatticeForceField, + AnalyticalFfromSuperF3D& analyticalExternalField, + T scale = T(1.) ); + ~ForceFromExtField3D() override {}; + void applyForce(typename std::deque >::iterator p, int pInt, + ParticleSystem3D& psSys) override; +private: +// SuperLattice3D& _sLattice; +// SuperLatticeField3D& _sLatticeForceField; + AnalyticalFfromSuperF3D& _analyticalExternalField; + T _scale; +}; + +} + +#endif diff --git a/src/particles/forces/forceFromExtField3D.hh b/src/particles/forces/forceFromExtField3D.hh new file mode 100755 index 0000000..1f9c567 --- /dev/null +++ b/src/particles/forces/forceFromExtField3D.hh @@ -0,0 +1,66 @@ +/* + * Copyright (C) 2018 Marie-Luise Maier + * E-mail contact: info@openlb.net + * The most recent release of OpenLB can be downloaded at + * + * + * 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. + */ + +#ifndef ForceFromExtField3D_HH +#define ForceFromExtField3D_HH + +#include "forceFromExtField3D.h" + +namespace olb { + +template class PARTICLETYPE, typename DESCRIPTOR> +ForceFromExtField3D::ForceFromExtField3D( +// SuperLattice3D& sLattice, + // SuperLatticeField3D& sLatticeForceField, + AnalyticalFfromSuperF3D& analyticalExternalField, + T scale) : + Force3D(), +// _sLattice(sLattice), + //_sLatticeForceField(sLatticeForceField), + _analyticalExternalField(analyticalExternalField), + _scale(scale) +{ } + +template class PARTICLETYPE, typename DESCRIPTOR> +void ForceFromExtField3D::applyForce( + typename std::deque >::iterator p, int pInt, + ParticleSystem3D& pSys) { + +// int glob = p->getCuboid(); +// int latticePos[3] = { 0, 0, 0 }; +// +// _sLattice.getCuboidGeometry().get(p->getCuboid()).getLatticeR(latticePos, +// &p->getPos()[0]); +// +// int X[4] = { glob, latticePos[0], latticePos[1], latticePos[2] }; +// + T F[3] = { T(), T(), T() }; +// _sLatticeForceField(F, X); + + _analyticalExternalField(F, &p->getPos()[0]); + p->getForce()[0] += F[0] * _scale; + p->getForce()[1] += F[1] * _scale; + p->getForce()[2] += F[2] * _scale; +} + +} +#endif diff --git a/src/particles/forces/forces.h b/src/particles/forces/forces.h new file mode 100644 index 0000000..e84d830 --- /dev/null +++ b/src/particles/forces/forces.h @@ -0,0 +1,42 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2016 Thomas Henn, Mathias J. Krause + * E-mail contact: info@openlb.net + * The most recent release of OpenLB can be downloaded at + * + * + * 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. +*/ + +/** \file + * Groups all the include .h-files for 3D particles forces in + * the particles forces directory. + */ + +#include "force3D.h" +#include "stokesDragForce3D.h" +#include "weightForce3D.h" +#include "buoyancyForce3D.h" +#include "hertzMindlinDeresiewicz3D.h" +#include "transferExternalForce3D.h" +#include "interpMagForceForMagP3D.h" +#include "magneticForceForMagP3D.h" +#include "linearDampingForceForMagDipoleMoment3D.h" +#include "stokesDragForceForHomVelField3D.h" +#include "magneticForceFromHField3D.h" +#include "forceFromExtField3D.h" + + diff --git a/src/particles/forces/forces.hh b/src/particles/forces/forces.hh new file mode 100644 index 0000000..52c0b90 --- /dev/null +++ b/src/particles/forces/forces.hh @@ -0,0 +1,40 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2016 Thomas Henn, Mathias J. Krause + * E-mail contact: info@openlb.net + * The most recent release of OpenLB can be downloaded at + * + * + * 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. +*/ + +/** \file + * Groups all the include .h-files for 3D particles forces in + * the particles forces directory. + */ + +#include "force3D.hh" +#include "stokesDragForce3D.hh" +#include "weightForce3D.hh" +#include "buoyancyForce3D.hh" +#include "hertzMindlinDeresiewicz3D.hh" +#include "transferExternalForce3D.hh" +#include "interpMagForceForMagP3D.hh" +#include "magneticForceForMagP3D.hh" +#include "linearDampingForceForMagDipoleMoment3D.hh" +#include "stokesDragForceForHomVelField3D.hh" +#include "magneticForceFromHField3D.hh" +#include "forceFromExtField3D.hh" diff --git a/src/particles/forces/hertzMindlinDeresiewicz3D.cpp b/src/particles/forces/hertzMindlinDeresiewicz3D.cpp new file mode 100644 index 0000000..f28c6b6 --- /dev/null +++ b/src/particles/forces/hertzMindlinDeresiewicz3D.cpp @@ -0,0 +1,36 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2016 Mathias J. Krause, Marie-Luise Maier, Sascha Janz + * E-mail contact: info@openlb.net + * The most recent release of OpenLB can be downloaded at + * + * + * 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. +*/ + +#include "dynamics/latticeDescriptors.h" + +#include "particles/particle3D.h" +#include "particles/particle3D.hh" +#include "hertzMindlinDeresiewicz3D.h" +#include "hertzMindlinDeresiewicz3D.hh" + +namespace olb { + +template class HertzMindlinDeresiewicz3D>; +template class HertzMindlinDeresiewicz3D>; + +} // namespace olb diff --git a/src/particles/forces/hertzMindlinDeresiewicz3D.h b/src/particles/forces/hertzMindlinDeresiewicz3D.h new file mode 100755 index 0000000..f839dd5 --- /dev/null +++ b/src/particles/forces/hertzMindlinDeresiewicz3D.h @@ -0,0 +1,75 @@ +/* + * 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 + * + * + * 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 + * + * validation paper for contact: + * H. Kruggel-Emden, E. Simsek, S. Rickelt, S. Wirtz, V. Scherer: Review and + * extension of normal force models for the Discrete Element Method, Powder + * Technology 171 (2007) 157-173 + */ + +#ifndef HERTZMINDLINDERESIEWICZ3D_H +#define HERTZMINDLINDERESIEWICZ3D_H + +#include +#include "functors/lattice/superLatticeLocalF3D.h" +#include "particles/particleSystem3D.h" +#include "force3D.h" + +namespace olb { + +template class PARTICLETYPE> +class ParticleSystem3D; + +template class PARTICLETYPE, typename DESCRIPTOR> +class HertzMindlinDeresiewicz3D: public Force3D { + +public: + HertzMindlinDeresiewicz3D(T G1, T G2, T v1, T v2, T scale1 = T(1.), T scale2 = T(1.), + bool validationKruggelEmden = false); + ~HertzMindlinDeresiewicz3D() override {}; + void applyForce(typename std::deque >::iterator p, int pInt, + ParticleSystem3D &pSys) override; + void computeForce(typename std::deque >::iterator p, int pInt, + ParticleSystem3D& pSys, T force[3]); +private: + T _G1; // Shear modulus (PA) + T _G2; // Shear modulus (PA) + T _v1; // poisson ratio + T _v2; // poisson ratio + T _scale1; // scales normal value of force + T _scale2; // scales tangential value of force + T E1, E2; // E-Modul Particle + T eE; // equivalent combined E-Modul + T eG; // equivalent combined E-Modul + // constant eta_n from paper H. Kruggel-Endem, Powder Technology 171 (2007) 157-173, + // Table 3 / brass particle + bool _validationKruggelEmden; // use values of paper +}; + +} + +#endif 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 + * + * + * 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 +#include "hertzMindlinDeresiewicz3D.h" + +namespace olb { + +template class PARTICLETYPE, typename DESCRIPTOR> +HertzMindlinDeresiewicz3D::HertzMindlinDeresiewicz3D( + T G1, T G2, T v1, T v2, T scale1, T scale2, bool validationKruggelEmden) : + Force3D(), _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 class PARTICLETYPE, typename DESCRIPTOR> +void HertzMindlinDeresiewicz3D::applyForce( + typename std::deque >::iterator p, int pInt, + ParticleSystem3D& pSys) +{ + T force[3] = {T(), T(), T()}; + + computeForce(p, pInt, pSys, force); +} + +template class PARTICLETYPE, typename DESCRIPTOR> +void HertzMindlinDeresiewicz3D::computeForce( + typename std::deque >::iterator p, int pInt, + ParticleSystem3D& pSys, T force[3]) +{ + + std::vector> ret_matches; + // kind of contactDetection has to be chosen in application + pSys.getContactDetection()->getMatches(pInt, ret_matches); + + PARTICLETYPE* 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 d_(_d); + Vector 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 diff --git a/src/particles/forces/hertzMindlinDeresiewiczForCombWithCollisionModel3D.cpp b/src/particles/forces/hertzMindlinDeresiewiczForCombWithCollisionModel3D.cpp new file mode 100644 index 0000000..f28c6b6 --- /dev/null +++ b/src/particles/forces/hertzMindlinDeresiewiczForCombWithCollisionModel3D.cpp @@ -0,0 +1,36 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2016 Mathias J. Krause, Marie-Luise Maier, Sascha Janz + * E-mail contact: info@openlb.net + * The most recent release of OpenLB can be downloaded at + * + * + * 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. +*/ + +#include "dynamics/latticeDescriptors.h" + +#include "particles/particle3D.h" +#include "particles/particle3D.hh" +#include "hertzMindlinDeresiewicz3D.h" +#include "hertzMindlinDeresiewicz3D.hh" + +namespace olb { + +template class HertzMindlinDeresiewicz3D>; +template class HertzMindlinDeresiewicz3D>; + +} // namespace olb diff --git a/src/particles/forces/hertzMindlinDeresiewiczForCombWithCollisionModel3D.h b/src/particles/forces/hertzMindlinDeresiewiczForCombWithCollisionModel3D.h new file mode 100755 index 0000000..d116c91 --- /dev/null +++ b/src/particles/forces/hertzMindlinDeresiewiczForCombWithCollisionModel3D.h @@ -0,0 +1,76 @@ +/* + * 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 + * + * + * 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 + * + * validation paper for contact: + * H. Kruggel-Emden, E. Simsek, S. Rickelt, S. Wirtz, V. Scherer: Review and + * extension of normal force models for the Discrete Element Method, Powder + * Technology 171 (2007) 157-173 + */ + +#ifndef HERTZMINDLINDERESIEWICZ3D_H +#define HERTZMINDLINDERESIEWICZ3D_H + +#include +#include "functors/lattice/superLatticeLocalF3D.h" +#include "particles/particleSystem3D.h" +#include "force3D.h" + +namespace olb { + +template class PARTICLETYPE> +class ParticleSystem3D; + +template class PARTICLETYPE, template< + typename W> class DESCRIPTOR> +class HertzMindlinDeresiewicz3D: public Force3D { + +public: + HertzMindlinDeresiewicz3D(T G1, T G2, T v1, T v2, T scale1 = T(1.), T scale2 = T(1.), + bool validationKruggelEmden = false); + ~HertzMindlinDeresiewicz3D() override {}; + void applyForce(typename std::deque >::iterator p, int pInt, + ParticleSystem3D &pSys) override ; + void computeForce(typename std::deque >::iterator p, int pInt, + ParticleSystem3D& pSys, T force[3]); +private: + T _G1; // Shear modulus (PA) + T _G2; // Shear modulus (PA) + T _v1; // poisson ratio + T _v2; // poisson ratio + T _scale1; // scales normal value of force + T _scale2; // scales tangential value of force + T E1, E2; // E-Modul Particle + T eE; // equivalent combined E-Modul + T eG; // equivalent combined E-Modul + // constant eta_n from paper H. Kruggel-Endem, Powder Technology 171 (2007) 157-173, + // Table 3 / brass particle + bool _validationKruggelEmden; // use values of paper +}; + +} + +#endif 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 + * + * + * 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 +#include "hertzMindlinDeresiewicz3D.h" + +namespace olb { + +template class PARTICLETYPE, template< + typename W> class DESCRIPTOR> +HertzMindlinDeresiewicz3D::HertzMindlinDeresiewicz3D( + T G1, T G2, T v1, T v2, T scale1, T scale2, bool validationKruggelEmden) : + Force3D(), _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 class PARTICLETYPE, template< + typename W> class DESCRIPTOR> +void HertzMindlinDeresiewicz3D::applyForce( + typename std::deque >::iterator p, int pInt, + ParticleSystem3D& pSys) +{ + T force[3] = {T(), T(), T()}; + + computeForce(p, pInt, pSys, force); +} + + +template class PARTICLETYPE, template< + typename W> class DESCRIPTOR> +void HertzMindlinDeresiewicz3D::computeForce( + typename std::deque >::iterator p, int pInt, + ParticleSystem3D& pSys, T force[3]) +{ + if (p->getSActivity() > 1) { + + std::vector> ret_matches; + // kind of contactDetection has to be chosen in application + pSys.getContactDetection()->getMatches(pInt, ret_matches); + + PARTICLETYPE* 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 d_(_d); + Vector 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 + + diff --git a/src/particles/forces/interpMagForceForMagP3D.cpp b/src/particles/forces/interpMagForceForMagP3D.cpp new file mode 100644 index 0000000..45a1c90 --- /dev/null +++ b/src/particles/forces/interpMagForceForMagP3D.cpp @@ -0,0 +1,40 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2015 Marie-Luise Maier, Mathias J. Krause + * E-mail contact: info@openlb.net + * The most recent release of OpenLB can be downloaded at + * + * + * 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. + */ + +/// Magnetic field that creates magnetization in wire has to be orthogonal to the wire. +/// to calculate the magnetic force on particles around a cylinder +/// (J. Lindner et al. / Computers and Chemical Engineering 54 (2013) 111-121) + +#include "dynamics/latticeDescriptors.h" + +#include "particles/particle3D.h" +#include "particles/particle3D.hh" +#include "interpMagForceForMagP3D.h" +#include "interpMagForceForMagP3D.hh" + +namespace olb { + +template class InterpMagForceForMagP3D>; + +} // namespace olb + diff --git a/src/particles/forces/interpMagForceForMagP3D.h b/src/particles/forces/interpMagForceForMagP3D.h new file mode 100644 index 0000000..4381509 --- /dev/null +++ b/src/particles/forces/interpMagForceForMagP3D.h @@ -0,0 +1,59 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2015 Marie-Luise Maier, Mathias J. Krause + * E-mail contact: info@openlb.net + * The most recent release of OpenLB can be downloaded at + * + * + * 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. + */ + +/// Magnetic field that creates magnetization in wire has to be orthogonal to the wire. +/// to calculate the magnetic force on particles around a cylinder +/// (J. Lindner et al. / Computers and Chemical Engineering 54 (2013) 111-121) +#ifndef InterpMagForceForMagP3D_H +#define InterpMagForceForMagP3D_H + +#include +#include "functors/lattice/superLatticeLocalF3D.h" +#include "particles/particleSystem3D.h" +#include "force3D.h" + +namespace olb { + +template class PARTICLETYPE> +class ParticleSystem3D; + +template class PARTICLETYPE, typename DESCRIPTOR> +class InterpMagForceForMagP3D : public Force3D { + +public: + InterpMagForceForMagP3D(T scale = T(1), T scaleT = T(1)); + InterpMagForceForMagP3D(InterpMagForceForMagP3D& f); + + ~InterpMagForceForMagP3D() override {}; + void applyForce(typename std::deque >::iterator p, + int pInt, ParticleSystem3D& psSys) override ; + +private: + T _scale; + T _scaleT; + T _sRad2; +}; + +} + +#endif // INTERPMAGFORCEFORMAGP3D_H_ diff --git a/src/particles/forces/interpMagForceForMagP3D.hh b/src/particles/forces/interpMagForceForMagP3D.hh new file mode 100644 index 0000000..a088bbb --- /dev/null +++ b/src/particles/forces/interpMagForceForMagP3D.hh @@ -0,0 +1,158 @@ +/* This file is part of the OpenLB library + * + * 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 + * + * + * 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. + */ + +/// Magnetic field that creates magnetization in wire has to be orthogonal to the wire. +/// to calculate the magnetic force on particles around a cylinder +/// (J. Lindner et al. / Computers and Chemical Engineering 54 (2013) 111-121) +#ifndef InterpMagForceForMagP3D_HH +#define InterpMagForceForMagP3D_HH + +#include +#include +#include "interpMagForceForMagP3D.h" + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +namespace olb { + +template class PARTICLETYPE, typename DESCRIPTOR> +InterpMagForceForMagP3D::InterpMagForceForMagP3D(T scale, T scaleT) + : Force3D(), + _scale(scale), + _scaleT(scaleT) +{ +// this->_name = "InterpMagForceForMagP3D"; +} + +template class PARTICLETYPE, typename DESCRIPTOR> +void InterpMagForceForMagP3D::applyForce( + typename std::deque >::iterator p, int pInt, + ParticleSystem3D& pSys) +{ + + std::vector < std::pair > ret_matches; + pSys.getDetection()->getMatches(pInt, ret_matches); + // std::random_shuffle ( ret_matches.begin(), ret_matches.end() ); + + /*const*/ PARTICLETYPE* p2 = NULL; + typename std::vector >::iterator it = + ret_matches.begin(); + + Vector dMom_1 = Vector((p->getMoment())); + if (dMom_1.norm() > std::numeric_limits::epsilon()) { + T m_p1 = p->getMagnetisation(); + + for (; it != ret_matches.end(); it++) { + if (it->second >= std::numeric_limits < T > ::epsilon()) { + + p2 = &pSys[it->first]; + + // No magnetic force between particles which are in physcal contact + // if ((p2->getRad() + p->getRad()) <= std::sqrt(it->second)) { + + // get neighbour particle moment + Vector dMom_2((p2->getMoment())); + if (dMom_2.norm() > std::numeric_limits::epsilon()) { + T m_p2 = p2->getMagnetisation(); + + // given moment magnitudes as scale factors + T m_i_scaleFactor = dMom_1.norm(); + T m_j_scaleFactor = dMom_2.norm(); + + // normalised moment directions + Vector n_i(dMom_1); + n_i.normalize(); + Vector n_j(dMom_2); + n_j.normalize(); + + // constants + T mu_0 = 4 * M_PI * 1.e-7; // magnetic constant + T mu_i = 4. / 3 * M_PI * pow(p->getRad(), 3) * m_p1 * m_i_scaleFactor ; // magnetic moment of particle i + T mu_j = 4. / 3 * M_PI * pow(p2->getRad(), 3) * m_p2 * m_j_scaleFactor ; // of particle j + + //vector from particle1 to particle2 + Vector r_ij; + //S Das wäre der Richtungsvektor von p2 zu p und nicht von p zu p2, aber die Rechung unten berücksichtigt das wohl, da richtige Ergebnisse + r_ij[0] = p->getPos()[0] - p2->getPos()[0]; + r_ij[1] = p->getPos()[1] - p2->getPos()[1]; + r_ij[2] = p->getPos()[2] - p2->getPos()[2]; + + // distance from particle1 to particle2 + T r = r_ij.norm(); + + // normalised direction vector + Vector t_ij(r_ij); + t_ij.normalize(); + + // FORCE of Dipole 1 on Dipole 2 + Vector mi = {n_i}; + mi *= mu_i; + Vector mj = {n_j}; + mj *= mu_j; + Vector rn = {r_ij}; + rn *= -1. ; + normalize(rn); + T scalar_termF = (3 * mu_0 ) / (4 * M_PI * std::pow(r, 4)); + Vector force = mj * (mi * rn) + mi * (mj * rn) + rn * (mi * mj) - 5 * rn * (mi * rn) * (mj * rn) ; + force *= scalar_termF; + + // TORQUE of Dipole 1 on Dipole 2 + // T_ij = -[(mu_0*mu_i*mu_j)/(4*M_PI*r^3)][n_i x n_j - 3(n_j * t_ij)n_i x t_ij] + T scalar_termT = - (mu_0 * mu_i * mu_j) / (4 * M_PI * std::pow(r, 3)); + Vector torque = crossProduct3D(n_i, n_j); + torque -= crossProduct3D( 3 * (n_j * t_ij) * n_i, t_ij); + torque *= scalar_termT; + + // Force an torque for overlapping particles + if ((p2->getRad() + p->getRad())