diff options
Diffstat (limited to 'src/particles/forces/interpMagForceForMagP3D.hh')
| -rw-r--r-- | src/particles/forces/interpMagForceForMagP3D.hh | 158 |
1 files changed, 158 insertions, 0 deletions
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 + * <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. + */ + +/// 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 <cmath> +#include <vector> +#include "interpMagForceForMagP3D.h" + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +namespace olb { + +template<typename T, template<typename U> class PARTICLETYPE, typename DESCRIPTOR> +InterpMagForceForMagP3D<T, PARTICLETYPE, DESCRIPTOR>::InterpMagForceForMagP3D(T scale, T scaleT) + : Force3D<T, PARTICLETYPE>(), + _scale(scale), + _scaleT(scaleT) +{ +// this->_name = "InterpMagForceForMagP3D"; +} + +template<typename T, template<typename U> class PARTICLETYPE, typename DESCRIPTOR> +void InterpMagForceForMagP3D<T, PARTICLETYPE, DESCRIPTOR>::applyForce( + typename std::deque<PARTICLETYPE<T> >::iterator p, int pInt, + ParticleSystem3D<T, PARTICLETYPE>& pSys) +{ + + std::vector < std::pair<size_t, T> > ret_matches; + pSys.getDetection()->getMatches(pInt, ret_matches); + // std::random_shuffle ( ret_matches.begin(), ret_matches.end() ); + + /*const*/ PARTICLETYPE<T>* p2 = NULL; + typename std::vector<std::pair<size_t, T> >::iterator it = + ret_matches.begin(); + + Vector<T, 3> dMom_1 = Vector<T, 3>((p->getMoment())); + if (dMom_1.norm() > std::numeric_limits<T>::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<T, 3> dMom_2((p2->getMoment())); + if (dMom_2.norm() > std::numeric_limits<T>::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<T, 3> n_i(dMom_1); + n_i.normalize(); + Vector<T, 3> 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<T, 3> 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, 3> t_ij(r_ij); + t_ij.normalize(); + + // FORCE of Dipole 1 on Dipole 2 + Vector<T, 3> mi = {n_i}; + mi *= mu_i; + Vector<T, 3> mj = {n_j}; + mj *= mu_j; + Vector<T, 3> rn = {r_ij}; + rn *= -1. ; + normalize(rn); + T scalar_termF = (3 * mu_0 ) / (4 * M_PI * std::pow(r, 4)); + Vector<T, 3> 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<T, 3> 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()) >= std::sqrt(it->second)) { + + normalize(force); + torque *= 0; + force *= mu_0 * pow((mu_i + mu_j) / 2., 2.) / (4 * M_PI * pow(r / 2, 4.)) ; + } + + p->getTorque()[0] += 0.5 * torque[0] * _scaleT; + p->getTorque()[1] += 0.5 * torque[1] * _scaleT; + p->getTorque()[2] += 0.5 * torque[2] * _scaleT; + p->getForce()[0] -= 0.5 * force[0] * _scale ; + p->getForce()[1] -= 0.5 * force[1] * _scale ; + p->getForce()[2] -= 0.5 * force[2] * _scale ; + p2->getTorque()[0] -= 0.5 * torque[0] * _scaleT; + p2->getTorque()[1] -= 0.5 * torque[1] * _scaleT; + p2->getTorque()[2] -= 0.5 * torque[2] * _scaleT; + p2->getForce()[0] += 0.5 * force[0] * _scale ; + p2->getForce()[1] += 0.5 * force[1] * _scale ; + p2->getForce()[2] += 0.5 * force[2] * _scale ; + + } + // } + } + } + } +} + +} + +#endif // INTERPMAGFORCEFORMAGP3D_H_ |