Initialize at openlb-1-3

1 files changed, 150 insertions, 0 deletions

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
+ *  <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.
+ */
+
+/* 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<typename T, typename Lattice, template<typename V> class Particle>
+DragModelBase<T,Lattice,Particle>::DragModelBase(UnitConverter<T, Lattice>& converter)
+         : _converter(converter)
+{}
+
+
+////////////////////// Class StokesSimplifiedDragModel ////////////////////////
+
+template<typename T, typename Lattice, template<typename V> class Particle>
+StokesSimplifiedDragModel<T,Lattice,Particle>::StokesSimplifiedDragModel(UnitConverter<T, Lattice>& converter)
+         : DragModelBase<T,Lattice,Particle>(converter)
+{}
+
+template<typename T, typename Lattice, template<typename V> class Particle>
+T StokesSimplifiedDragModel<T,Lattice,Particle>::operator() (
+                         Particle<T>* p, T latticeVelF[], T latticeVelP[], int globicFull[] )
+{
+  return 1.83;
+}
+
+
+////////////////////// Class MorsiDragModel ////////////////////////
+
+template<typename T, typename Lattice, template<typename V> class Particle>
+MorsiDragModel<T,Lattice,Particle>::MorsiDragModel(UnitConverter<T, Lattice>& converter)
+         : DragModelBase<T,Lattice,Particle>(converter)
+{
+  this->_reP = std::make_shared<NewtonianParticleReynoldsNumber<T,Lattice,Particle> > (this->_converter);
+}
+
+template<typename T, typename Lattice, template<typename V> class Particle>
+T MorsiDragModel<T,Lattice,Particle>::operator() (
+                         Particle<T>* 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<typename T, typename Lattice, template<typename V> class Particle>
+DewsburyDragModel<T,Lattice,Particle>::DewsburyDragModel(UnitConverter<T, Lattice>& converter)
+         : DragModelBase<T,Lattice,Particle>(converter)
+{
+  this->_reP = std::make_shared<NewtonianParticleReynoldsNumber<T,Lattice,Particle> > (this->_converter);
+}
+
+template<typename T, typename Lattice, template<typename V> class Particle>
+T DewsburyDragModel<T,Lattice,Particle>::operator() (
+                         Particle<T>* 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<typename T, typename Lattice, template<typename V> class Particle>
+PowerLawDewsburyDragModel<T,Lattice,Particle>::PowerLawDewsburyDragModel (
+        UnitConverter<T, Lattice>& converter, SuperLattice3D<T, Lattice>& sLattice )
+         : DewsburyDragModel<T,Lattice,Particle>(converter)
+{
+  this->_reP = std::make_shared<PowerLawParticleReynoldsNumber<T,Lattice,Particle> > (this->_converter, sLattice);
+}
+
+
+}
+
+#endif