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/twoWayCouplings/backCouplingModels.h | 124 ++++++ .../twoWayCouplings/backCouplingModels.hh | 230 +++++++++++ src/particles/twoWayCouplings/dragModels3D.h | 113 ++++++ src/particles/twoWayCouplings/dragModels3D.hh | 150 +++++++ .../twoWayCouplings/forwardCouplingModels3D.h | 148 +++++++ .../twoWayCouplings/forwardCouplingModels3D.hh | 280 +++++++++++++ src/particles/twoWayCouplings/twoWayCouplings3D.h | 32 ++ src/particles/twoWayCouplings/twoWayCouplings3D.hh | 32 ++ .../twoWayCouplings/twoWayHelperFunctionals.h | 232 +++++++++++ .../twoWayCouplings/twoWayHelperFunctionals.hh | 431 +++++++++++++++++++++ 10 files changed, 1772 insertions(+) create mode 100644 src/particles/twoWayCouplings/backCouplingModels.h create mode 100644 src/particles/twoWayCouplings/backCouplingModels.hh create mode 100644 src/particles/twoWayCouplings/dragModels3D.h create mode 100644 src/particles/twoWayCouplings/dragModels3D.hh create mode 100644 src/particles/twoWayCouplings/forwardCouplingModels3D.h create mode 100644 src/particles/twoWayCouplings/forwardCouplingModels3D.hh create mode 100644 src/particles/twoWayCouplings/twoWayCouplings3D.h create mode 100644 src/particles/twoWayCouplings/twoWayCouplings3D.hh create mode 100644 src/particles/twoWayCouplings/twoWayHelperFunctionals.h create mode 100644 src/particles/twoWayCouplings/twoWayHelperFunctionals.hh (limited to 'src/particles/twoWayCouplings') diff --git a/src/particles/twoWayCouplings/backCouplingModels.h b/src/particles/twoWayCouplings/backCouplingModels.h new file mode 100644 index 0000000..8f016d2 --- /dev/null +++ b/src/particles/twoWayCouplings/backCouplingModels.h @@ -0,0 +1,124 @@ +/* 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 back-coupling methods -- header file. + */ + +#ifndef LB_BACK_COUPLING_MODELS_H +#define LB_BACK_COUPLING_MODELS_H + +#include "twoWayHelperFunctionals.h" + +namespace olb { + +/** Abstact base class for BaseBackCouplingModel. + * Its raison d'etre consists of not being templetized in Lattice. + */ +template class Particle> +class BackCouplingModel { +public: + /// Class operator to apply the coupling, for overload. + virtual bool operator() (Particle* p, int globic, int material, int subCycles=1)=0; + /// Resets external field + virtual void resetExternalField(int material)=0; +}; + +/** Abstact class for all the back-coupling models, + * viz., momentum coupling from particle to fluid. + * Input parameters in attice units. + */ +template class Particle> +class BaseBackCouplingModel : public BackCouplingModel { +public: + /// Resets external field + virtual void resetExternalField(int material) override; +protected: + /// Constructor + BaseBackCouplingModel ( UnitConverter& converter, + SuperLattice3D& sLattice, + SuperGeometry3D& sGeometry ); + UnitConverter& _converter; // reference to a UnitConverter + SuperGeometry3D& _sGeometry; + SuperLattice3D& _sLattice; // reference to a lattice +private: + // Pointers to functions to reset fluid force + std::shared_ptr > _zeroAnalytical; + std::shared_ptr > _zeroField; +}; + +/** Back-coupling is performed on the cell containing the particle + * and its neighbours within a cube of one lattice spacing, as per in Maier et al. (2017). + * Input parameters in attice units. + */ +template class Particle> +class CubicDeltaBackCouplingModel : public BaseBackCouplingModel { +public: + /// Constructor + CubicDeltaBackCouplingModel ( UnitConverter& converter, + SuperLattice3D& sLattice, + SuperGeometry3D& sGeometry ); + /// Class operator to apply the coupling. + virtual bool operator() (Particle* p, int globic, int material, int subCycles=1) override; +protected: + int _range = 1; + T _delta[4][4][4] = { T() }; + std::shared_ptr > _cubicDeltaFunctional; +}; + +/** Back-coupling is performed only on the cell containing the particle. + * Input parameters in attice units. + */ +template class Particle> +class LocalBackCouplingModel : public BaseBackCouplingModel { +public: + /// Constructor + LocalBackCouplingModel ( UnitConverter& converter, + SuperLattice3D& sLattice, + SuperGeometry3D& sGeometry ); + /// Class operator to apply the coupling. + virtual bool operator() (Particle* p, int globic, int material, int subCycles=1) override; +}; + +/** Class for a generic non-local back-coupling model, + * viz., momentum coupling from particle to fluid. + * It reproduces the characteristics (viz., smoothing) of an input forward coupling model. + * Input parameters in attice units. + */ +template class Particle> +class NonLocalBaseBackCouplingModel : public BaseBackCouplingModel { +public: + /// Constructor + NonLocalBaseBackCouplingModel ( UnitConverter& converter, + SuperLattice3D& sLattice, + SuperGeometry3D& sGeometry, + SmoothingFunctional& smoothingFunctional ); + /// Class operator to apply the coupling. + virtual bool operator() (Particle* p, int globic, int material, int subCycles=1) override; +protected: + SmoothingFunctional& _smoothingFunctional; // Functional to treat non-local smoothing +}; + +} + +#endif diff --git a/src/particles/twoWayCouplings/backCouplingModels.hh b/src/particles/twoWayCouplings/backCouplingModels.hh new file mode 100644 index 0000000..e88c62a --- /dev/null +++ b/src/particles/twoWayCouplings/backCouplingModels.hh @@ -0,0 +1,230 @@ +/* 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 back-coupling methods -- generic implementation. + */ + +#ifndef LB_BACK_COUPLING_MODELS_HH +#define LB_BACK_COUPLING_MODELS_HH + +namespace olb { + + +////////////////////// Class BaseBackCouplingModel //////////////////////// + +template class Particle> +BaseBackCouplingModel::BaseBackCouplingModel ( + UnitConverter& converter, + SuperLattice3D& sLattice, + SuperGeometry3D& sGeometry ) + : _converter(converter), + _sGeometry(sGeometry), + _sLattice(sLattice) +{ + _zeroAnalytical = std::make_shared > (T()); + _zeroField = std::make_shared > (*_zeroAnalytical, *_zeroAnalytical, *_zeroAnalytical); +} + +template class Particle> +void BaseBackCouplingModel::resetExternalField(int material) +{ + // resets external field + this->_sLattice.template defineField(this->_sGeometry, material, *_zeroField); + + // NECESSARY to communicate values before using them in operator() + this->_sLattice.communicate(); +} + + +////////////////////// Class CubicDeltaBackCouplingModel //////////////////////// + +template class Particle> +CubicDeltaBackCouplingModel::CubicDeltaBackCouplingModel ( + UnitConverter& converter, + SuperLattice3D& sLattice, + SuperGeometry3D& sGeometry ) + : BaseBackCouplingModel(converter, sLattice, sGeometry) +{ + _cubicDeltaFunctional = std::make_shared > ( + this->_sLattice, this->_converter, this->_sGeometry ); +} + +template class Particle> +bool CubicDeltaBackCouplingModel::operator() (Particle* p, int globic, int material, int subCycles) +{ + int locIC = this->_sLattice.getLoadBalancer().loc(globic); + + // reading the force from the value stored inside the particle + std::vector physForceP = p->getStoreForce(); // physical force acting on the particle + + T latticeForceP[3] = {T(), T(), T()}; // dimensionless force acting on the particle + latticeForceP[0] = physForceP[0] / this->_converter.getConversionFactorForce(); + latticeForceP[1] = physForceP[1] / this->_converter.getConversionFactorForce(); + latticeForceP[2] = physForceP[2] / this->_converter.getConversionFactorForce(); + + 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 ); + + // smooth Dirac delta + this->_cubicDeltaFunctional->operator() (_delta, physPosP, globic); + + T tempDelta = T(); + T F[3] = {T(), T(), T()}; // dimensionless smoothed force + + for (int i = -_range; i <= _range; ++i) { + for (int j = -_range; j <= _range; ++j) { + for (int k = -_range; k <= _range; ++k) { + if (this->_sGeometry.getBlockGeometry(locIC).getMaterial( + latticeRoundedPosP[0] + i, latticeRoundedPosP[1] + j, + latticeRoundedPosP[2] + k) == material) { + + tempDelta = _delta[i + _range][j + _range][k + _range]; + + F[0] = -latticeForceP[0] * tempDelta / (T)(subCycles); + F[1] = -latticeForceP[1] * tempDelta / (T)(subCycles); + F[2] = -latticeForceP[2] * tempDelta / (T)(subCycles); + + this->_sLattice.getBlockLattice(locIC).get ( + latticeRoundedPosP[0] + i, + latticeRoundedPosP[1] + j, + latticeRoundedPosP[2] + k ).template addField( F ); + } + } + } + } + return true; +} + + +////////////////////// Class LocalDeltaBackCouplingModel //////////////////////// + +template class Particle> +LocalBackCouplingModel::LocalBackCouplingModel ( + UnitConverter& converter, + SuperLattice3D& sLattice, + SuperGeometry3D& sGeometry ) + : BaseBackCouplingModel(converter, sLattice, sGeometry) +{} + +template class Particle> +bool LocalBackCouplingModel::operator() (Particle* p, int globic, int material, int subCycles) +{ + int locIC = this->_sLattice.getLoadBalancer().loc(globic); + + // reading the force from the value stored inside the particle + std::vector physForceP = p->getStoreForce(); // physical force acting on the particle + + T latticeForceP[3] = {T(), T(), T()}; // dimensionless force acting on the particle + latticeForceP[0] = physForceP[0] / this->_converter.getConversionFactorForce(); + latticeForceP[1] = physForceP[1] / this->_converter.getConversionFactorForce(); + latticeForceP[2] = physForceP[2] / this->_converter.getConversionFactorForce(); + + 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 ); + + if (this->_sGeometry.getBlockGeometry(locIC).getMaterial( + latticeRoundedPosP[0], latticeRoundedPosP[1], + latticeRoundedPosP[2]) == material) { + + T F[3] = {T(), T(), T()}; // dimensionless smoothed force + F[0] = -latticeForceP[0] / (T)(subCycles); + F[1] = -latticeForceP[1] / (T)(subCycles); + F[2] = -latticeForceP[2] / (T)(subCycles); + + this->_sLattice.getBlockLattice(locIC).get ( + latticeRoundedPosP[0], + latticeRoundedPosP[1], + latticeRoundedPosP[2] ).template addField( F ); + } + + return true; +} + + +////////////////////// Class NonLocalBaseBackCouplingModel //////////////////////// + +template class Particle> +NonLocalBaseBackCouplingModel::NonLocalBaseBackCouplingModel ( + UnitConverter& converter, + SuperLattice3D& sLattice, + SuperGeometry3D& sGeometry, + SmoothingFunctional& smoothingFunctional ) + : BaseBackCouplingModel(converter, sLattice, sGeometry), + _smoothingFunctional(smoothingFunctional) +{} + +template class Particle> +bool NonLocalBaseBackCouplingModel::operator() (Particle* p, int globic, int material, int subCycles) +{ + int locIC = this->_sLattice.getLoadBalancer().loc(globic); + + // reading the force from the value stored inside the particle + std::vector physForceP = p->getStoreForce(); // physical force acting on the particle + T latticeForceP[3] = {T(), T(), T()}; // dimensionless force acting on the particle + latticeForceP[0] = physForceP[0] / this->_converter.getConversionFactorForce(); + latticeForceP[1] = physForceP[1] / this->_converter.getConversionFactorForce(); + latticeForceP[2] = physForceP[2] / this->_converter.getConversionFactorForce(); + + // Updating force through voxels within kernel smoothing length from the bubble's position + for (int i=0; i_smoothingFunctional.getSize(); i++) { + + // Position of the iterated voxel + int iLatticePosF[3] = {0, 0, 0}; + this->_smoothingFunctional.getLatticePos(iLatticePosF, i); + + // Updating iterated voxel + if (this->_sGeometry.getBlockGeometry(locIC).getMaterial( + iLatticePosF[0], iLatticePosF[1], + iLatticePosF[2]) == material) { + + // Weighted force acting on the iterated voxel + T F[3] = {T(), T(), T()}; // dimensionless smoothed force + F[0] = -latticeForceP[0] * this->_smoothingFunctional.getWeight(i) / (T)(subCycles); + F[1] = -latticeForceP[1] * this->_smoothingFunctional.getWeight(i) / (T)(subCycles); + F[2] = -latticeForceP[2] * this->_smoothingFunctional.getWeight(i) / (T)(subCycles); + + this->_sLattice.getBlockLattice(locIC).get ( + iLatticePosF[0], iLatticePosF[1], iLatticePosF[2] ).template addField( F ); + } + } + return true; +} + +} + +#endif diff --git a/src/particles/twoWayCouplings/dragModels3D.h b/src/particles/twoWayCouplings/dragModels3D.h new file mode 100644 index 0000000..0f22429 --- /dev/null +++ b/src/particles/twoWayCouplings/dragModels3D.h @@ -0,0 +1,113 @@ +/* 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. + */ + +/* Drag force models for Lagrangian two-way coupling -- header file. + */ + +#ifndef LB_DRAG_MODELS_H +#define LB_DRAG_MODELS_H + +#include "functors/lattice/reductionF3D.h" +#include "twoWayHelperFunctionals.h" + +namespace olb { + +/** Abstact base class for DragModelBase. + * Its raison d'etre consists of not being templetized in Lattice. + */ +template class Particle> +class DragModel { +public: + /// Returns the scalar drag coefficient to overload. globicFull = { globic, latticeRoundedP[0, ..., 2] } + virtual T operator() (Particle* p, T latticeVelF[], T latticeVelP[], int globicFull[])=0; +protected: + /// Functional to compute particle Reynolds number + std::shared_ptr > _reP; +}; + +/** Abstact class for all the drag models. + * The virtual method computeDragCoeff returns the drag coefficient. + * Input parameters in attice units. + */ +template class Particle> +class DragModelBase : public DragModel { +public: + /// Constructor + DragModelBase(UnitConverter& converter); + /// Returns the scalar drag coefficient to overload. + //virtual T operator() (Particle* p, T latticeVelF[], T latticeVelP[], int globic)=0; +protected: + UnitConverter& _converter; // reference to a UnitConverter +}; + +/** Class to compute a drag coefficient Cd=1.83 for low-Re Stokes drag. + */ +template class Particle> +class StokesSimplifiedDragModel : public DragModelBase { +public: + /// Constructor + StokesSimplifiedDragModel(UnitConverter& converter); + /// Returns the scalar drag coefficient. globicFull = { globic, latticeRoundedP[0, ..., 2] } + virtual T operator() (Particle* p, T latticeVelF[], T latticeVelP[], int globicFull[]) override; +}; + +/** Class to compute the standard drag coefficient + * as in Morsi and Alexander (1972). + */ +template class Particle> +class MorsiDragModel : public DragModelBase { +public: + /// Constructor + MorsiDragModel(UnitConverter& converter); + /// Returns the scalar drag coefficient. globicFull = { globic, latticeRoundedP[0, ..., 2] } + virtual T operator() (Particle* p, T latticeVelF[], T latticeVelP[], int globicFull[]) override; +}; + +/** Class to compute the drag coefficient for gas bubbles in a liquid fluid phase + * as in Dewsbury et al. (1999). + */ +template class Particle> +class DewsburyDragModel : public DragModelBase { +public: + /// Constructor + DewsburyDragModel(UnitConverter& converter); + /// Returns the scalar drag coefficient. globicFull = { globic, latticeRoundedP[0, ..., 2] } + virtual T operator() (Particle* p, T latticeVelF[], T latticeVelP[], int globicFull[]) override; +}; + +/** Class to compute the drag coefficient for gas bubbles in a liquid fluid phase + * as in Dewsbury et al. (1999), in a power-law fluid. + */ +template class Particle> +class PowerLawDewsburyDragModel : public DewsburyDragModel { +public: + /// Constructor + PowerLawDewsburyDragModel ( + UnitConverter& converter, SuperLattice3D& sLattice ); +}; + + +} + +#endif diff --git a/src/particles/twoWayCouplings/dragModels3D.hh b/src/particles/twoWayCouplings/dragModels3D.hh new file mode 100644 index 0000000..0027544 --- /dev/null +++ b/src/particles/twoWayCouplings/dragModels3D.hh @@ -0,0 +1,150 @@ +/* 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. + */ + +/* Drag force models for Lagrangian two-way coupling -- generic implementation. + */ + +#ifndef LB_DRAG_MODELS_HH +#define LB_DRAG_MODELS_HH + +namespace olb { + + +////////////////////// Class DragModelBase //////////////////////// + +template class Particle> +DragModelBase::DragModelBase(UnitConverter& converter) + : _converter(converter) +{} + + +////////////////////// Class StokesSimplifiedDragModel //////////////////////// + +template class Particle> +StokesSimplifiedDragModel::StokesSimplifiedDragModel(UnitConverter& converter) + : DragModelBase(converter) +{} + +template class Particle> +T StokesSimplifiedDragModel::operator() ( + Particle* p, T latticeVelF[], T latticeVelP[], int globicFull[] ) +{ + return 1.83; +} + + +////////////////////// Class MorsiDragModel //////////////////////// + +template class Particle> +MorsiDragModel::MorsiDragModel(UnitConverter& converter) + : DragModelBase(converter) +{ + this->_reP = std::make_shared > (this->_converter); +} + +template class Particle> +T MorsiDragModel::operator() ( + Particle* p, T latticeVelF[], T latticeVelP[], int globicFull[] ) +{ + T physVelRelative = this->_converter.getPhysVelocity ( + sqrt( pow(latticeVelF[0] - latticeVelP[0],2) + + pow(latticeVelF[1] - latticeVelP[1],2) + + pow(latticeVelF[2] - latticeVelP[2],2) ) ); + + T ReP = this->_reP->operator() (p, physVelRelative, globicFull); + + T a[3] = {T(), T(), T()}; + if (ReP < 0.1) { + a[0] = 0.0; a[1] = 24.0; a[2] = 0.0; + } + else if (ReP < 1.0) { + a[0] = 3.69; a[1] = 22.73; a[2] = 0.0903; + } + else if (ReP < 10.0) { + a[0] = 1.222; a[1] = 29.16667; a[2] =-3.8889; + } + else if (ReP < 100.0) { + a[0] = 0.6167; a[1] = 46.5; a[2] =-116.67; + } + else if (ReP < 1000.0) { + a[0] = 0.3644; a[1] = 498.33; a[2] =-2778; + } + else if (ReP < 5000.0) { + a[0] = 0.357; a[1] = 148.62; a[2] =-4.75e4; + } + else if (ReP < 10000.0) { + a[0] = 0.46; a[1] =-490.546; a[2] = 57.87e4; + } + else { + a[0] = 0.5191; a[1] =-1662.5; a[2] = 5.4167e6; + } + + + return ( a[0] + a[1]/ReP + a[2]/(ReP*ReP) ) * physVelRelative; +} + + +////////////////////// Class DewsburyDragModel //////////////////////// + +template class Particle> +DewsburyDragModel::DewsburyDragModel(UnitConverter& converter) + : DragModelBase(converter) +{ + this->_reP = std::make_shared > (this->_converter); +} + +template class Particle> +T DewsburyDragModel::operator() ( + Particle* p, T latticeVelF[], T latticeVelP[], int globicFull[] ) +{ + T physVelRelative = this->_converter.getPhysVelocity ( + sqrt( pow(latticeVelF[0] - latticeVelP[0],2) + + pow(latticeVelF[1] - latticeVelP[1],2) + + pow(latticeVelF[2] - latticeVelP[2],2) ) ); + + T ReP = this->_reP->operator() (p, physVelRelative, globicFull); + + T Cd = 0.95; + if (ReP <= 195.) { + Cd = 16./ReP * (1. + 0.173*pow(ReP, 0.657)) + + 0.413 / (1. + 16300*pow(ReP, -1.09)); + } + return Cd * this->_converter.getLatticeVelocity(physVelRelative); +} + + +////////////////////// Class PowerLawDewsburyDragModel //////////////////////// + +template class Particle> +PowerLawDewsburyDragModel::PowerLawDewsburyDragModel ( + UnitConverter& converter, SuperLattice3D& sLattice ) + : DewsburyDragModel(converter) +{ + this->_reP = std::make_shared > (this->_converter, sLattice); +} + + +} + +#endif diff --git a/src/particles/twoWayCouplings/forwardCouplingModels3D.h b/src/particles/twoWayCouplings/forwardCouplingModels3D.h new file mode 100644 index 0000000..17ff63c --- /dev/null +++ b/src/particles/twoWayCouplings/forwardCouplingModels3D.h @@ -0,0 +1,148 @@ +/* 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 -- header file. + */ + +#ifndef LB_FORWARD_COUPLING_MODELS_H +#define LB_FORWARD_COUPLING_MODELS_H + +#include "functors/lattice/reductionF3D.h" +#include "twoWayHelperFunctionals.h" + +namespace olb { + +/** Abstact base class for LocalBaseCouplingModel. + * Its raison d'etre consists of not being templetized in Lattice. + */ +template class Particle> +class ForwardCouplingModel { +public: + /// Class operator to apply the coupling, for overload. + virtual bool operator() (Particle* p, int globic)=0; +}; + +/** Abstact class for all the local forward-coupling models, + * viz., momentum coupling from fluid to particle. + * Input parameters in attice units. + */ +template class Particle> +class LocalBaseForwardCouplingModel : public ForwardCouplingModel { +public: + /// Class operator to apply the coupling, for overload. + virtual bool operator() (Particle* p, int globic) override; +protected: + /// Constructor + LocalBaseForwardCouplingModel ( UnitConverter& converter, + SuperLattice3D& sLattice, + SuperGeometry3D& sGeometry, + DragModel& dragModel ); + SuperGeometry3D& _sGeometry; + UnitConverter& _converter; // reference to a UnitConverter + SuperLattice3D& _sLattice; // reference to a lattice + DragModel& _dragModel; // reference to a drag model + // Functional to interpolate lattice density at particle's location + std::shared_ptr > _interpLatticeDensity; + // Functional to interpolate lattice velocity at particle's location + std::shared_ptr > _interpLatticeVelocity; + // Momentum-exchange helper functional + std::shared_ptr > _momentumExchange; +}; + +/** Class for a naive forward-coupling model. + * Input parameters in attice units. + */ +template class Particle> +class NaiveForwardCouplingModel : public LocalBaseForwardCouplingModel { +public: + /// Constructor + NaiveForwardCouplingModel ( UnitConverter& converter, + SuperLattice3D& sLattice, + SuperGeometry3D& sGeometry, + DragModel& dragModel ); +}; + +/** Class for a forward-coupling model following the Ladd's mechanism. + * Input parameters in attice units. + */ +template class Particle> +class LaddForwardCouplingModel : public LocalBaseForwardCouplingModel { +public: + /// Constructor + LaddForwardCouplingModel ( UnitConverter& converter, + SuperLattice3D& sLattice, + SuperGeometry3D& sGeometry, + DragModel& dragModel ); +}; + +/** Abstact class for all the non-local forward-coupling models, + * viz., momentum coupling from fluid to particle. + * Input parameters in attice units. + */ +template class Particle> +class NonLocalBaseForwardCouplingModel : public ForwardCouplingModel { +public: + /// Class operator to apply the coupling, for overload. + virtual bool operator() (Particle* p, int globic) override; +protected: + /// Constructor + NonLocalBaseForwardCouplingModel ( UnitConverter& converter, + SuperLattice3D& sLattice, + SuperGeometry3D& sGeometry, + DragModel& dragModel, + SmoothingFunctional& smoothingFunctional ); + SuperGeometry3D& _sGeometry; + UnitConverter& _converter; // reference to a UnitConverter + SuperLattice3D& _sLattice; // reference to a lattice + DragModel& _dragModel; // reference to a drag model + SmoothingFunctional& _smoothingFunctional; // Functional to treat non-local smoothing + // Functional to interpolate lattice density at particle's location. + // It is assumed that it returns the exact value for the exact cell's location! + std::shared_ptr > _interpLatticeDensity; + // Functional to interpolate lattice velocity at particle's location + // It is assumed that it returns the exact value for the exact cell's location! + std::shared_ptr > _interpLatticeVelocity; + // Momentum-exchange helper functional + std::shared_ptr > _momentumExchange; +}; + +/** Class for a naive, non-local forward-coupling model as in Sungkorn et al. (2011), but with + * an extra-normalization of the smoothing function. + * Input parameters in attice units. + */ +template class Particle> +class NaiveNonLocalForwardCouplingModel : public NonLocalBaseForwardCouplingModel { +public: + /// Constructor + NaiveNonLocalForwardCouplingModel ( UnitConverter& converter, + SuperLattice3D& sLattice, + SuperGeometry3D& sGeometry, + DragModel& dragModel, + SmoothingFunctional& smoothingFunctional ); +}; + + +} + +#endif 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 diff --git a/src/particles/twoWayCouplings/twoWayCouplings3D.h b/src/particles/twoWayCouplings/twoWayCouplings3D.h new file mode 100644 index 0000000..c0bdce8 --- /dev/null +++ b/src/particles/twoWayCouplings/twoWayCouplings3D.h @@ -0,0 +1,32 @@ +/* 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. + */ + +/* \file + * Groups all the 3D include files in the twoWayCoupling folder. + */ + +#include "twoWayHelperFunctionals.h" +#include "dragModels3D.h" +#include "forwardCouplingModels3D.h" +#include "backCouplingModels.h" diff --git a/src/particles/twoWayCouplings/twoWayCouplings3D.hh b/src/particles/twoWayCouplings/twoWayCouplings3D.hh new file mode 100644 index 0000000..1d8cdc1 --- /dev/null +++ b/src/particles/twoWayCouplings/twoWayCouplings3D.hh @@ -0,0 +1,32 @@ +/* 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. + */ + +/* \file + * Groups all the generic 3D template files in the twoWayCoupling folder. + */ + +#include "twoWayHelperFunctionals.hh" +#include "dragModels3D.hh" +#include "forwardCouplingModels3D.hh" +#include "backCouplingModels.hh" diff --git a/src/particles/twoWayCouplings/twoWayHelperFunctionals.h b/src/particles/twoWayCouplings/twoWayHelperFunctionals.h new file mode 100644 index 0000000..7e76138 --- /dev/null +++ b/src/particles/twoWayCouplings/twoWayHelperFunctionals.h @@ -0,0 +1,232 @@ +/* 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. + */ + +/* Helper functionals for Lagrangian two-way coupling methods -- header file. + */ + +#ifndef LB_TWO_WAY_HELPER_FUNCTIONALS_H +#define LB_TWO_WAY_HELPER_FUNCTIONALS_H + +#include "functors/lattice/reductionF3D.h" + +namespace olb { + +/** Data structure for smoothing functionals. + * Stores the lattice position of a cell within smoothing kernel length + * and the related multiplicative weight. + */ +template +struct LatticePosAndWeight { + int globic = 0; + int latticePos[3] = {0, 0, 0}; + T weight = T(); +}; + +/** Abstract class for particle Reynolds number computation within drag model. + * Its raison d'etre consists of not being templetized in Lattice. + */ +template class Particle> +class ParticleReynoldsNumber { +public: + /// Returns the particle Reynolds number. globicFull = { globic, latticeRoundedP[0, ..., 2] } + virtual T operator() (Particle* p, T magU, int globicFull[])=0; + /// Destructor + virtual ~ParticleReynoldsNumber() {}; +protected: + T _RePmin = 0.01; +}; + +/** Abstract class for particle Reynolds number computation within drag model. + */ +template class Particle> +class ParticleReynoldsNumberBase : public ParticleReynoldsNumber { +public: + /// Destructor + virtual ~ParticleReynoldsNumberBase() {}; +protected: + /// Constructor + ParticleReynoldsNumberBase(UnitConverter& converter); + UnitConverter& _converter; // reference to a UnitConverter +}; + +/** Class class for Newtonian particle Reynolds number computation within drag model. + */ +template class Particle> +class NewtonianParticleReynoldsNumber: public ParticleReynoldsNumberBase { +public: + /// Constructor + NewtonianParticleReynoldsNumber(UnitConverter& converter); + /// Destructor + ~NewtonianParticleReynoldsNumber() {}; + /// Returns the particle Reynolds number. globicFull = { globic, latticeRoundedP[0, ..., 2] } + virtual T operator() (Particle* p, T magU, int globicFull[]) override; +}; + +/** Class class for power-law particle Reynolds number computation within drag model. + */ +template class Particle> +class PowerLawParticleReynoldsNumber: public ParticleReynoldsNumberBase { +public: + /// Constructor + PowerLawParticleReynoldsNumber ( + UnitConverter& converter, SuperLattice3D& sLattice ); + /// Destructor + ~PowerLawParticleReynoldsNumber() {}; + /// Returns the particle Reynolds number. globicFull = { globic, latticeRoundedP[0, ..., 2] } + virtual T operator() (Particle* p, T magU, int globicFull[]) override; +protected: + SuperLattice3D& _sLattice; // reference to a lattice +}; + +/** Abstact class for all the local forward-coupling models, + * viz., momentum coupling from fluid to particle. + * Input parameters in attice units. + */ +template +class TwoWayHelperFunctional { +public: + /// Computes the momentum transfer from fluid to particle. + virtual bool operator() ( T gF[], T latticeVelF[], T latticeVelP[], + T physPosP[], int latticeRoundedP[], + int globic )=0; + virtual ~TwoWayHelperFunctional(); +protected: + /// Constructor + TwoWayHelperFunctional ( UnitConverter& converter, + SuperLattice3D& sLattice ); + UnitConverter& _converter; // reference to a UnitConverter + SuperLattice3D& _sLattice; // reference to a lattice + std::shared_ptr > _interpLatticeDensity; + std::shared_ptr > _interpLatticeVelocity ; +}; + +/** Naive way + */ +template +class NaiveMomentumExchange : public TwoWayHelperFunctional { +public: + /// Constructor + NaiveMomentumExchange ( UnitConverter& converter, + SuperLattice3D& sLattice, + std::shared_ptr > interpLatticeDensity ); + /// Computes the momentum transfer from fluid to particle. + virtual bool operator() ( T gF[], T latticeVelF[], T latticeVelP[], + T physPosP[], int latticeRoundedP[], + int globic ) override; +}; + +/** Using Ladd mechanism + */ +template +class LaddMomentumExchange : public TwoWayHelperFunctional { +public: + /// Constructor + LaddMomentumExchange ( UnitConverter& converter, + SuperLattice3D& sLattice, + std::shared_ptr > interpLatticeDensity, + std::shared_ptr > interpLatticeVelocity ); + /// Computes the momentum transfer from fluid to particle. + virtual bool operator() ( T gF[], T latticeVelF[], T latticeVelP[], + T physPosP[], int latticeRoundedP[], + int globic ) override; +}; + +/** Abstact class for all the smoothing functionals. + */ +template +class SmoothingFunctional { +public: + // Returns the size of _latticePosAndWeight + int getSize(); + // Returns the lattice position of the i-th element of _latticePosAndWeight + void getLatticePos(int latticePos[], int i); + // Returns the globic of the i-th element of _latticePosAndWeight + int getGlobic(int i); + // Returns the weight relative to the i-th element of _latticePosAndWeight + T getWeight(int i); + // Rebuilds _latticePosAndWeight with the new cells within _kernelLength from the bubble's position + bool update(T physPosP[], int globic); +protected: + /// Constructor + SmoothingFunctional(T kernelLength, UnitConverter& converter, SuperLattice3D& sLattice); + /// The actual smoothing function + virtual T smoothingFunction(T delta)=0; + /// Returns the weight for smoothing. + virtual T compute(T physPosP[], T physPosL[])=0; + T _kernelLength; // Kernel's smoothing length. + UnitConverter& _converter; // reference to a UnitConverter + SuperLattice3D& _sLattice; // reference to a lattice + // positions and weights of the cells within _kernelLength from bubble's position + std::deque > _latticePosAndWeight; +}; + +/** Abstact class for all the linear-averaging smoothing functionals. + */ +template +class LinearAveragingSmoothingFunctional : public SmoothingFunctional { +protected: + /// Constructor + LinearAveragingSmoothingFunctional(T kernelLength, UnitConverter& converter, SuperLattice3D& sLattice); + /// Returns the weight for smoothing. + virtual T compute(T physPosP[], T physPosL[]) override; +}; + +/** Abstact class for all the volume-averaging smoothing functionals. + */ +template +class VolumeAveragingSmoothingFunctional : public SmoothingFunctional { +protected: + /// Constructor + VolumeAveragingSmoothingFunctional(T kernelLength, UnitConverter& converter, SuperLattice3D& sLattice); + /// Returns the weight for smoothing. + virtual T compute(T physPosP[], T physPosL[]) override; +}; + +/** Smoothing functional as in Deen et al (2004), Chem. Eng. Sci 59. + */ +template +class DeenSmoothingFunctional : public LinearAveragingSmoothingFunctional { +public: + /// Constructor + DeenSmoothingFunctional(T kernelLength, UnitConverter& converter, SuperLattice3D& sLattice); +protected: + /// The actual smoothing function + virtual T smoothingFunction(T delta) override; +}; + +/** Stepwise smoothing functional. + */ +template +class StepSmoothingFunctional : public VolumeAveragingSmoothingFunctional { +public: + /// Constructor + StepSmoothingFunctional(T kernelLength, UnitConverter& converter, SuperLattice3D& sLattice); +protected: + /// The actual smoothing function + virtual T smoothingFunction(T delta) override; +}; + +} + +#endif diff --git a/src/particles/twoWayCouplings/twoWayHelperFunctionals.hh b/src/particles/twoWayCouplings/twoWayHelperFunctionals.hh new file mode 100644 index 0000000..b7db1b3 --- /dev/null +++ b/src/particles/twoWayCouplings/twoWayHelperFunctionals.hh @@ -0,0 +1,431 @@ +/* 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. + */ + +/* Helper functionals for Lagrangian two-way coupling methods -- generic implementation. + */ + +#ifndef LB_TWO_WAY_HELPER_FUNCTIONALS_HH +#define LB_TWO_WAY_HELPER_FUNCTIONALS_HH + +namespace olb { + +////////////////////// Class ParticleReynoldsNumberBase //////////////////////// + +template class Particle> +ParticleReynoldsNumberBase::ParticleReynoldsNumberBase(UnitConverter& converter) + : _converter(converter) +{} + + +////////////////////// Class NewtonianParticleReynoldsNumber //////////////////////// + +template class Particle> +NewtonianParticleReynoldsNumber::NewtonianParticleReynoldsNumber(UnitConverter& converter) + : ParticleReynoldsNumberBase(converter) +{} + +template class Particle> +T NewtonianParticleReynoldsNumber::operator() ( Particle* p, T magU, int globicFull[]) +{ + T ReP = 2. * p->getRad() * magU / this->_converter.getPhysViscosity(); + return ReP > this->_RePmin ? ReP : this->_RePmin; +} + + +////////////////////// Class PowerLawParticleReynoldsNumber //////////////////////// + +template class Particle> +PowerLawParticleReynoldsNumber::PowerLawParticleReynoldsNumber ( + UnitConverter& converter, SuperLattice3D& sLattice ) + : ParticleReynoldsNumberBase(converter), + _sLattice(sLattice) +{} + +template class Particle> +T PowerLawParticleReynoldsNumber::operator() ( Particle* p, T magU, int globicFull[]) +{ + // loc() indicates the local cuboid number locIC of the actual processing thread, + // for given global cuboid number iC + // this is to get appropriate particle system on locIC + int locIC = _sLattice.getLoadBalancer().loc(globicFull[0]); + + T physPosP[3] = {T(), T(), T()}; // particle's physical position + physPosP[0