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 --- .../twoWayCouplings/forwardCouplingModels3D.hh | 280 +++++++++++++++++++++ 1 file changed, 280 insertions(+) create mode 100644 src/particles/twoWayCouplings/forwardCouplingModels3D.hh (limited to 'src/particles/twoWayCouplings/forwardCouplingModels3D.hh') diff --git a/src/particles/twoWayCouplings/forwardCouplingModels3D.hh b/src/particles/twoWayCouplings/forwardCouplingModels3D.hh new file mode 100644 index 0000000..b07d50b --- /dev/null +++ b/src/particles/twoWayCouplings/forwardCouplingModels3D.hh @@ -0,0 +1,280 @@ +/* Lattice Boltzmann sample, written in C++, using the OpenLB + * library + * + * Copyright (C) 2019 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. + */ + +/* Models for Lagrangian forward-coupling methods -- generic implementation. + */ + +#ifndef LB_FORWARD_COUPLING_MODELS_HH +#define LB_FORWARD_COUPLING_MODELS_HH + +namespace olb { + + +////////////////////// Class LocalBaseCouplingModel //////////////////////// + +template class Particle> +LocalBaseForwardCouplingModel::LocalBaseForwardCouplingModel ( + UnitConverter& converter, + SuperLattice3D& sLattice, + SuperGeometry3D& sGeometry, + DragModel& dragModel ) + : _sGeometry(sGeometry), + _converter(converter), + _sLattice(sLattice), + _dragModel(dragModel) +{ + this->_interpLatticeDensity = std::make_shared > ( + this->_sLattice, _sGeometry, this->_converter ); + this->_interpLatticeVelocity = std::make_shared > ( + this->_sLattice, this->_converter ); +} + +template class Particle> +bool LocalBaseForwardCouplingModel::operator() (Particle* p, int globic) +{ + /// Getting the particle and its containing cell's position + T physPosP[3] = {T(), T(), T()}; // particle's physical position + physPosP[0] = (p->getPos()[0]); + physPosP[1] = (p->getPos()[1]); + physPosP[2] = (p->getPos()[2]); + // particle's dimensionless position, rounded at neighbouring voxel + int latticeRoundedPosP[3] = {0, 0, 0}; + this->_sLattice.getCuboidGeometry().get(globic).getLatticeR ( + latticeRoundedPosP, physPosP ); + // { globic, latticeRoundedP[0, ..., 2] } + int globicFull[4] = { globic, + latticeRoundedPosP[0], + latticeRoundedPosP[1], + latticeRoundedPosP[2] }; + + // Particle's velocity + T physVelP[3] = {T(), T(), T()}; // Physical + physVelP[0] = (p->getVel()[0]); + physVelP[1] = (p->getVel()[1]); + physVelP[2] = (p->getVel()[2]); + // Lattice + T latticeVelP[3] = {T(), T(), T()}; // particle's dimensionless velocity + latticeVelP[0] = this->_converter.getLatticeVelocity(physVelP[0]); + latticeVelP[1] = this->_converter.getLatticeVelocity(physVelP[1]); + latticeVelP[2] = this->_converter.getLatticeVelocity(physVelP[2]); + + // Lattice's velocity at particle's location + T physVelF[3] = {T(), T(), T()}; // Physical + this->_interpLatticeVelocity->operator() (physVelF, physPosP, globic); + // Lattice + T latticeVelF[3] = {T(), T(), T()}; // Lattice's dimensionless velocity at particle's location + latticeVelF[0] = this->_converter.getLatticeVelocity(physVelF[0]); + latticeVelF[1] = this->_converter.getLatticeVelocity(physVelF[1]); + latticeVelF[2] = this->_converter.getLatticeVelocity(physVelF[2]); + + // Computing fluid-particle momentum transfer + T gF[3] = {T(), T(), T()}; // force density gF + this->_momentumExchange->operator() ( gF, latticeVelF, latticeVelP, + physPosP, latticeRoundedPosP, globic); + + // Computing drag coefficient + T Cd = this->_dragModel(p, latticeVelF, latticeVelP, globicFull); + + /// Computing drag force in dimensionless units + T latticePRad = p->getRad() / _converter.getConversionFactorLength(); + T latticeForceP[3] = {T(), T(), T()}; // dimensionless force acting on the particle + latticeForceP[0] = .5 * Cd * M_PI*pow(latticePRad,2) * gF[0] * (latticeVelF[0] - latticeVelP[0]); + latticeForceP[1] = .5 * Cd * M_PI*pow(latticePRad,2) * gF[1] * (latticeVelF[1] - latticeVelP[1]); + latticeForceP[2] = .5 * Cd * M_PI*pow(latticePRad,2) * gF[2] * (latticeVelF[2] - latticeVelP[2]); + + /// Computing physical drag force + std::vector physForceP(3, T()); // physical force acting on the particle + physForceP[0] = latticeForceP[0] * this->_converter.getConversionFactorForce(); + physForceP[1] = latticeForceP[1] * this->_converter.getConversionFactorForce(); + physForceP[2] = latticeForceP[2] * this->_converter.getConversionFactorForce(); + + /// Updating the particle + p->setStoreForce(physForceP); + return true; +} + + +////////////////////// Class NaiveForwardCouplingModel //////////////////////// + +template class Particle> +NaiveForwardCouplingModel::NaiveForwardCouplingModel ( + UnitConverter& converter, + SuperLattice3D& sLattice, + SuperGeometry3D& sGeometry, + DragModel& dragModel ) + : LocalBaseForwardCouplingModel(converter, sLattice, sGeometry, dragModel) +{ + this->_momentumExchange = std::make_shared > ( + this->_converter, this->_sLattice, this->_interpLatticeDensity ); +} + + +////////////////////// Class LaddForwardCouplingModel //////////////////////// + +template class Particle> +LaddForwardCouplingModel::LaddForwardCouplingModel ( + UnitConverter& converter, + SuperLattice3D& sLattice, + SuperGeometry3D& sGeometry, + DragModel& dragModel ) + : LocalBaseForwardCouplingModel(converter, sLattice, sGeometry, dragModel) +{ + this->_momentumExchange = std::make_shared > ( + this->_converter, this->_sLattice, + this->_interpLatticeDensity, this->_interpLatticeVelocity ); +} + + +////////////////////// Class NonLocalBaseCouplingModel //////////////////////// + +template class Particle> +NonLocalBaseForwardCouplingModel::NonLocalBaseForwardCouplingModel ( + UnitConverter& converter, + SuperLattice3D& sLattice, + SuperGeometry3D& sGeometry, + DragModel& dragModel, + SmoothingFunctional& smoothingFunctional ) + : _sGeometry(sGeometry), + _converter(converter), + _sLattice(sLattice), + _dragModel(dragModel), + _smoothingFunctional(smoothingFunctional) +{ + this->_interpLatticeDensity = std::make_shared > ( + this->_sLattice, _sGeometry, this->_converter ); + this->_interpLatticeVelocity = std::make_shared > ( + this->_sLattice, this->_converter ); +} + +template class Particle> +bool NonLocalBaseForwardCouplingModel::operator() (Particle* p, int globic) +{ + /// Getting the particle and its containing cell's position + T physPosP[3] = {T(), T(), T()}; // particle's physical position + physPosP[0] = (p->getPos()[0]); + physPosP[1] = (p->getPos()[1]); + physPosP[2] = (p->getPos()[2]); + //std::cout << "globic=" << globic << " physPosP=(" << physPosP[0] << ", " << physPosP[1] << ", " << physPosP[2] << ")" << std::endl; // <--- + // particle's dimensionless position, rounded at neighbouring voxel + int latticeRoundedPosP[3] = {0, 0, 0}; + this->_sLattice.getCuboidGeometry().get(globic).getLatticeR ( + latticeRoundedPosP, physPosP ); + // { globic, latticeRoundedP[0, ..., 2] } + int globicFull[4] = { globic, + latticeRoundedPosP[0], + latticeRoundedPosP[1], + latticeRoundedPosP[2] }; + + // Particle's velocity + T physVelP[3] = {T(), T(), T()}; // Physical + physVelP[0] = (p->getVel()[0]); + physVelP[1] = (p->getVel()[1]); + physVelP[2] = (p->getVel()[2]); + // Lattice + T latticeVelP[3] = {T(), T(), T()}; // particle's dimensionless velocity + latticeVelP[0] = this->_converter.getLatticeVelocity(physVelP[0]); + latticeVelP[1] = this->_converter.getLatticeVelocity(physVelP[1]); + latticeVelP[2] = this->_converter.getLatticeVelocity(physVelP[2]); + //std::cout << "globic=" << globic << " latticeVelP=(" << latticeVelP[0] << ", " << latticeVelP[1] << ", " << latticeVelP[2] << ")" << std::endl; // <--- + + // Update the smoothing functional + if ( ! this->_smoothingFunctional.update(physPosP, globic)) { + std::cout << "ERROR: no lattice point enclosed in particle's kernel length!" << std::endl; + return false; + } + + /// Computing dimensionless drag force through smoothed average within the kernel smoothing length + T latticeForceP[3] = {T(), T(), T()}; // dimensionless force acting on the particle + for (int i=0; i_smoothingFunctional.getSize(); i++) { + + // Position of the iterated voxel + int iLatticePosF[3] = {0, 0, 0}; + this->_smoothingFunctional.getLatticePos(iLatticePosF, i); + // Physical + T iPhysPosF[3] = {T(), T(), T()}; + this->_sLattice.getCuboidGeometry().get(globic).getPhysR ( + iPhysPosF, iLatticePosF ); + + // Fluid velocity at the iterated voxel + T iPhysVelF[3] = {T(), T(), T()}; // Physical + this->_interpLatticeVelocity->operator() (iPhysVelF, iPhysPosF, globic); + //std::cout << "globic=" << globic << " iPhysVelF=(" << iPhysVelF[0] << ", " << iPhysVelF[1] << ", " << iPhysVelF[2] << ")" << std::endl; // <--- + // Lattice + T iLatticeVelF[3] = {T(), T(), T()}; // Lattice's dimensionless velocity at particle's location + iLatticeVelF[0] = this->_converter.getLatticeVelocity(iPhysVelF[0]); + iLatticeVelF[1] = this->_converter.getLatticeVelocity(iPhysVelF[1]); + iLatticeVelF[2] = this->_converter.getLatticeVelocity(iPhysVelF[2]); + + // Computing fluid-particle momentum transfer + T gF[3] = {T(), T(), T()}; // force density gF + this->_momentumExchange->operator() ( gF, iLatticeVelF, latticeVelP, + physPosP, iLatticePosF, globic); + //std::cout << "globic=" << globic << " gF=(" << gF[0] << ", " << gF[1] << ", " << gF[2] << ")" << std::endl; // <--- + + // Computing drag coefficient + T Cd = this->_dragModel(p, iLatticeVelF, latticeVelP, globicFull); + //std::cout << "globic=" << globic << " Cd=" << Cd << std::endl; //<--- + + /// Computing drag force in dimensionless units + T latticePRad = p->getRad() / _converter.getConversionFactorLength(); + latticeForceP[0] += .5 * Cd * M_PI*pow(latticePRad,2) * gF[0] * (iLatticeVelF[0] - latticeVelP[0]) + * this->_smoothingFunctional.getWeight(i); + latticeForceP[1] += .5 * Cd * M_PI*pow(latticePRad,2) * gF[1] * (iLatticeVelF[1] - latticeVelP[1]) + * this->_smoothingFunctional.getWeight(i); + latticeForceP[2] += .5 * Cd * M_PI*pow(latticePRad,2) * gF[2] * (iLatticeVelF[2] - latticeVelP[2]) + * this->_smoothingFunctional.getWeight(i); + } + + /// Computing physical drag force + std::vector physForceP(3, T()); // physical force acting on the particle + physForceP[0] = latticeForceP[0] * this->_converter.getConversionFactorForce(); + physForceP[1] = latticeForceP[1] * this->_converter.getConversionFactorForce(); + physForceP[2] = latticeForceP[2] * this->_converter.getConversionFactorForce(); + + /// Updating the particle + p->setStoreForce(physForceP); + //std::cout << "globic=" << globic << " physForceP=(" << physForceP[0] << ", " << physForceP[1] << ", " << physForceP[2] << ")" << std::endl; // <--- + return true; +} + + +////////////////////// Class NaiveNonLocalForwardCouplingModel //////////////////////// + +template class Particle> +NaiveNonLocalForwardCouplingModel::NaiveNonLocalForwardCouplingModel ( + UnitConverter& converter, + SuperLattice3D& sLattice, + SuperGeometry3D& sGeometry, + DragModel& dragModel, + SmoothingFunctional& smoothingFunctional ) + : NonLocalBaseForwardCouplingModel(converter, sLattice, sGeometry, dragModel, smoothingFunctional) +{ + this->_momentumExchange = std::make_shared > ( + this->_converter, this->_sLattice, this->_interpLatticeDensity ); +} + + +} + +#endif -- cgit v1.2.3