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/hlbm/hlbmDynamics2D.hh | 261 +++++++++++++++++++++++++++++++++++
1 file changed, 261 insertions(+)
create mode 100644 src/particles/hlbm/hlbmDynamics2D.hh
(limited to 'src/particles/hlbm/hlbmDynamics2D.hh')
diff --git a/src/particles/hlbm/hlbmDynamics2D.hh b/src/particles/hlbm/hlbmDynamics2D.hh
new file mode 100644
index 0000000..95e996f
--- /dev/null
+++ b/src/particles/hlbm/hlbmDynamics2D.hh
@@ -0,0 +1,261 @@
+/* DESCRIPTOR Boltzmann sample, written in C++, using the OpenLB
+ * library
+ *
+ * Copyright (C) 2006-2016 Thomas Henn, Fabian Klemens, Robin Trunk, Davide Dapelo
+ * 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 PARTICLEDYNAMICS_2D_HH
+#define PARTICLEDYNAMICS_2D_HH
+
+#include "functors/analytical/indicator/indicCalc2D.h"
+#include "functors/lattice/superLatticeLocalF2D.h"
+#include "hlbmDynamics2D.h"
+
+namespace olb {
+
+template
+void ParticleDynamics2D::addCircle(Vector< T, 2> center, T radius, T density, T epsilon, Vector vel)
+{
+ _vectorOfIndicator.push_back (
+ new SmoothIndicatorCircle2D(center, radius, epsilon, density, vel) );
+}
+
+template
+void ParticleDynamics2D::addCuboid(Vector< T, 2> center, T xLength, T yLength, T density, T epsilon, T theta, Vector vel)
+{
+ _vectorOfIndicator.push_back (
+ new SmoothIndicatorCuboid2D(center, xLength, yLength, epsilon, theta, density, vel) );
+}
+
+template
+void ParticleDynamics2D::addTriangle(Vector< T, 2> center, T radius, T density, T epsilon, T theta, Vector vel)
+{
+ _vectorOfIndicator.push_back (
+ new SmoothIndicatorTriangle2D(center, radius, epsilon, theta, density, vel) );
+}
+
+template
+void ParticleDynamics2D::addParticle(SmoothIndicatorF2D& indicator)
+{
+ _vectorOfIndicator.push_back(&indicator);
+}
+
+template
+void ParticleDynamics2D::computeBoundaryForce(std::vector* >& indicator)
+{
+ SuperLatticePorousMomentumLossForce2D force(_sLattice, _superGeometry, indicator,_converter);
+ T sumF[force.getTargetDim()];
+ for (int i=0; i* >::size_type iInd=0; iInd!=indicator.size(); iInd++) {
+ /// get particle acceleration through boundary force and gravity (and buoyancy)
+ Vector acceleration2;
+ Vector force;
+ force[0] = sumF[0+4*iInd];
+ force[1] = sumF[1+4*iInd];
+ acceleration2[0] = sumF[0+4*iInd] / indicator[iInd]->getMass() + _accExt[0];
+ acceleration2[1] = sumF[1+4*iInd] / indicator[iInd]->getMass() + _accExt[1];
+ indicator[iInd]->setForce( force );
+ indicator[iInd]->setAcc2( acceleration2 );
+ indicator[iInd]->setAlpha2( sumF[2+4*iInd] / indicator[iInd]->getMofi() );
+ }
+}
+
+template
+void ParticleDynamics2D::addWallColl(SmoothIndicatorF2D& indicator, T delta)
+{
+ T w1 = 1e-7 / 2.;
+ T w = 1e-7 / 2.;
+
+ T rad = indicator.getRadius();
+ T massInv = 1. / indicator.getMass();
+
+ std::vector dx(2, T());
+ dx[0] = _lengthX - indicator.getPos()[0];
+ dx[1] = _lengthY - indicator.getPos()[1];
+
+ for (int i = 0; i < 2; i++) {
+ if (dx[i] <= rad) {
+ indicator.getAcc2()[i] += massInv * -dx[i] * (rad - dx[i]) / w1;
+ }
+ if (indicator.getPos()[i] <= rad) {
+ indicator.getAcc2()[i] += massInv * indicator.getPos()[i] * (rad - indicator.getPos()[i]) / w1;
+ }
+ if (dx[i] > rad && dx[i] <= rad + delta) {
+ indicator.getAcc2()[i] += massInv * -dx[i] * std::pow((rad + delta - dx[i]), 2) / w;
+ }
+ if (indicator.getPos()[i] > rad && indicator.getPos()[i] <= rad + delta) {
+ indicator.getAcc2()[i] += massInv * indicator.getPos()[i] * std::pow((rad + delta - indicator.getPos()[i]), 2) / w;
+ }
+ }
+}
+
+template
+void ParticleDynamics2D::verletIntegration(SmoothIndicatorF2D& indicator)
+{
+ T time = _converter.getConversionFactorTime();
+ T time2 = time*time;
+
+ Vector position, velocity;
+ Vector rotationMatrix;
+ for (int i=0; i<2; i++) {
+ position[i] = indicator.getPos()[i] + indicator.getVel()[i] * time + (0.5 * indicator.getAcc()[i] * time2);
+ T avgAcc = (indicator.getAcc()[i] + indicator.getAcc2()[i]) * 0.5;
+ velocity[i] = indicator.getVel()[i] + avgAcc * time;
+ }
+ indicator.setPos(position);
+ indicator.setVel(velocity);
+ indicator.setAcc(indicator.getAcc2());
+
+ indicator.setTheta( std::fmod( indicator.getTheta() + indicator.getOmega() * time + (0.5 * indicator.getAlpha() * time2), 2.*M_PI ) );
+ T avgAlpha = (indicator.getAlpha() + indicator.getAlpha2()) * 0.5;
+ indicator.setOmega( indicator.getOmega() + avgAlpha * time);
+ indicator.setAlpha( indicator.getAlpha2() );
+
+ T cos = std::cos(indicator.getTheta());
+ T sin = std::sin(indicator.getTheta());
+
+ rotationMatrix[0] = cos;
+ rotationMatrix[1] = -sin;
+ rotationMatrix[2] = sin;
+ rotationMatrix[3] = cos;
+ indicator.setRotationMatrix(rotationMatrix);
+}
+
+template
+void ParticleDynamics2D::updateParticleDynamics(std::string name, SmoothIndicatorF2D& indicator)
+{
+ this->verletIntegration(indicator);
+}
+
+template
+void ParticleDynamics2D::checkAndRemoveEscaped()
+{
+ for (auto i=_vectorOfIndicator.begin(); i!=_vectorOfIndicator.end(); i++) {
+ auto internal = std::make_shared > (*_indicatorF.get(), -_converter.getConversionFactorLength()+(**i).getEpsilon() );
+ IndicMinus2D layer(_indicatorF, internal);
+ SuperLatticeIndicatorSmoothIndicatorIntersection2D intersection(
+ _sLattice, _superGeometry, layer, **i );
+ int input[1];
+ T output[1] = {0.};
+ intersection(output, input);
+ if (output[0] == 1) {
+ _vectorOfIndicator.erase(i);
+ if (i==_vectorOfIndicator.end() ) {
+ break;
+ }
+ }
+ }
+}
+
+
+template
+void ParticleDynamics2D::addParticleField(SmoothIndicatorF2D& indicator)
+{
+ /// Analytical2D functor for particle motion (trans+rot)
+ ParticleU2D velocity(indicator, _converter);
+ _sLattice.setExternalParticleField(_superGeometry, velocity, indicator);
+}
+
+template
+void ParticleDynamics2D::simulateTimestep(std::string name)
+{
+ // Remove particles from domain if _escapeFromDomain is toggled
+ if (_escapeFromDomain) checkAndRemoveEscaped();
+
+ // Compute force acting on particles boundary
+ computeBoundaryForce(_vectorOfIndicator);
+
+ // Update particle dynamics and porous particle field
+ for (auto i=_vectorOfIndicator.begin(); i!=_vectorOfIndicator.end(); i++) {
+ updateParticleDynamics(name, **i);
+ addParticleField(**i);
+ }
+}
+
+template
+void ParticleDynamics2D::print()
+{
+ OstreamManager clout(std::cout, "ParticleDynamics2D");
+ clout << "Number of particles=" << _vectorOfIndicator.size() << std::endl;
+ int count = 0;
+ for (auto i=_vectorOfIndicator.begin(); i!=_vectorOfIndicator.end(); i++) {
+ clout << "Particle " << ++count << " (" << (**i).name() << "):" << std::endl;
+ clout << " |Circum radius(m)= " << std::setw(13) << (**i).getCircumRadius() <> name;
+
+ if (name == "circle") {
+ T radius = T(), mass = T();
+ Vector center = {T(), T()};
+ Vector vel = {T(), T()};
+ iout >> center[0] >> center[1] >> radius >> mass >> vel[0] >> vel[1];
+ addCircle(center, radius, mass, epsilon, vel);
+ }
+ }
+
+ iout.close();
+}
+
+// WORKS ONLY FOR CIRCLES FOR NOW!!
+template
+void ParticleDynamics2D::save(std::string filename)
+{
+ std::ofstream fout( (singleton::directories().getLogOutDir() + filename).c_str() );
+
+ for (auto i=_vectorOfIndicator.begin(); i!=_vectorOfIndicator.end(); i++) {
+ if ( (**i).name() == "circle") {
+ fout << (**i).name() << " "
+ << (**i).getPos()[0] << " " << (**i).getPos()[1] << " "
+ << (**i).getCircumRadius() << " " << (**i).getMass() << " "
+ << (**i).getVel()[0] << " " << (**i).getVel()[1]
+ << std::endl;
+ }
+ }
+
+ fout.close();
+}
+
+}
+
+#endif
--
cgit v1.2.3