diff options
Initialize at openlb-1-3
Diffstat (limited to 'src/particles/hlbm/hlbmDynamics2D.hh')
-rw-r--r-- | src/particles/hlbm/hlbmDynamics2D.hh | 261 |
1 files changed, 261 insertions, 0 deletions
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 + * <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. + */ + + + +#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<typename T, typename DESCRIPTOR> +void ParticleDynamics2D<T, DESCRIPTOR>::addCircle(Vector< T, 2> center, T radius, T density, T epsilon, Vector<S,2> vel) +{ + _vectorOfIndicator.push_back ( + new SmoothIndicatorCircle2D<T,T,true>(center, radius, epsilon, density, vel) ); +} + +template<typename T, typename DESCRIPTOR> +void ParticleDynamics2D<T, DESCRIPTOR>::addCuboid(Vector< T, 2> center, T xLength, T yLength, T density, T epsilon, T theta, Vector<S,2> vel) +{ + _vectorOfIndicator.push_back ( + new SmoothIndicatorCuboid2D<T,T,true>(center, xLength, yLength, epsilon, theta, density, vel) ); +} + +template<typename T, typename DESCRIPTOR> +void ParticleDynamics2D<T, DESCRIPTOR>::addTriangle(Vector< T, 2> center, T radius, T density, T epsilon, T theta, Vector<S,2> vel) +{ + _vectorOfIndicator.push_back ( + new SmoothIndicatorTriangle2D<T,T,true>(center, radius, epsilon, theta, density, vel) ); +} + +template<typename T, typename DESCRIPTOR> +void ParticleDynamics2D<T, DESCRIPTOR>::addParticle(SmoothIndicatorF2D<T,T,true>& indicator) +{ + _vectorOfIndicator.push_back(&indicator); +} + +template<typename T, typename DESCRIPTOR> +void ParticleDynamics2D<T, DESCRIPTOR>::computeBoundaryForce(std::vector<SmoothIndicatorF2D<T,T,true>* >& indicator) +{ + SuperLatticePorousMomentumLossForce2D<T, DESCRIPTOR> force(_sLattice, _superGeometry, indicator,_converter); + T sumF[force.getTargetDim()]; + for (int i=0; i<force.getTargetDim(); i++) { + sumF[i]=0.; + } + int input[1]; + force(sumF, input); + for (typename std::vector<SmoothIndicatorF2D<T,T,true>* >::size_type iInd=0; iInd!=indicator.size(); iInd++) { + /// get particle acceleration through boundary force and gravity (and buoyancy) + Vector<T,2> acceleration2; + Vector<T,2> 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<typename T, typename DESCRIPTOR> +void ParticleDynamics2D<T, DESCRIPTOR>::addWallColl(SmoothIndicatorF2D<T,T,true>& indicator, T delta) +{ + T w1 = 1e-7 / 2.; + T w = 1e-7 / 2.; + + T rad = indicator.getRadius(); + T massInv = 1. / indicator.getMass(); + + std::vector<T> 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<typename T, typename DESCRIPTOR> +void ParticleDynamics2D<T, DESCRIPTOR>::verletIntegration(SmoothIndicatorF2D<T,T,true>& indicator) +{ + T time = _converter.getConversionFactorTime(); + T time2 = time*time; + + Vector<T,2> position, velocity; + Vector<T,4> 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<typename T, typename DESCRIPTOR> +void ParticleDynamics2D<T, DESCRIPTOR>::updateParticleDynamics(std::string name, SmoothIndicatorF2D<T,T,true>& indicator) +{ + this->verletIntegration(indicator); +} + +template<typename T, typename DESCRIPTOR> +void ParticleDynamics2D<T, DESCRIPTOR>::checkAndRemoveEscaped() +{ + for (auto i=_vectorOfIndicator.begin(); i!=_vectorOfIndicator.end(); i++) { + auto internal = std::make_shared<IndicatorLayer2D<T> > (*_indicatorF.get(), -_converter.getConversionFactorLength()+(**i).getEpsilon() ); + IndicMinus2D<T> layer(_indicatorF, internal); + SuperLatticeIndicatorSmoothIndicatorIntersection2D<T,DESCRIPTOR,true> 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<typename T, typename DESCRIPTOR> +void ParticleDynamics2D<T, DESCRIPTOR>::addParticleField(SmoothIndicatorF2D<T,T,true>& indicator) +{ + /// Analytical2D functor for particle motion (trans+rot) + ParticleU2D<T,T,DESCRIPTOR> velocity(indicator, _converter); + _sLattice.setExternalParticleField(_superGeometry, velocity, indicator); +} + +template<typename T, typename DESCRIPTOR> +void ParticleDynamics2D<T, DESCRIPTOR>::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<typename T, typename DESCRIPTOR> +void ParticleDynamics2D<T, DESCRIPTOR>::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() <<std::endl; + clout << " |Mass(kg)= " << std::setw(13) << (**i).getMass() << std::endl; + clout << " |Position(m)= (" << std::setw(13) << (**i).getPos()[0] << ", " << std::setw(13) << (**i).getPos()[1] << ")" << std::endl; + clout << " |Angle(°)= " << std::setw(13) << (**i).getTheta()*(180/M_PI) << std::endl; + clout << " |Velocity(m/s)= (" << std::setw(13) <<(**i).getVel()[0] << ", " << std::setw(13) << (**i).getVel()[1] << ")" << std::endl; + clout << " |Ang. velocity(°/s)= " << std::setw(13) << (**i).getOmega()*(180/M_PI) << std::endl; + clout << " |Hydro. Force(N)= (" << std::setw(13) <<(**i).getForce()[0] << ", " << std::setw(13) << (**i).getForce()[1] << ")" << std::endl; + clout << " |Acceleration(m/s^2)= (" << std::setw(13) << (**i).getAcc()[0] << ", " << std::setw(13) << (**i).getAcc()[1] << ")" << std::endl; + clout << " |Ang. acc.(°/s^2)= " << std::setw(13) << (**i).getAlpha()*(180/M_PI) << std::endl; + } +} + +template<typename T, typename DESCRIPTOR> +void ParticleDynamics2D<T, DESCRIPTOR>::load(std::string filename, T epsilon) +{ + std::ifstream iout( (singleton::directories().getLogOutDir() + filename).c_str() ); + + while (iout) { + std::string name = ""; + iout >> name; + + if (name == "circle") { + T radius = T(), mass = T(); + Vector<T, 2> center = {T(), T()}; + Vector<T, 2> 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<typename T, typename DESCRIPTOR> +void ParticleDynamics2D<T, DESCRIPTOR>::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 |