Added functionality, tests and showcase for discrete particle method coupling

apps/showcases/BidisperseFluidizedBed/CMakeLists.txt

+waLBerla_add_executable ( NAME BidisperseFluidizedBedDPM
+                          DEPENDS blockforest boundary core domain_decomposition field lbm pe pe_coupling postprocessing timeloop vtk )
\ No newline at end of file

apps/showcases/CMakeLists.txt

+add_subdirectory( BidisperseFluidizedBed )
\ No newline at end of file

src/pe_coupling/all.h

@@ -24,6 +24,7 @@
 
 #pragma once
 
+#include "discrete_particle_methods/all.h"
 #include "mapping/all.h"
 #include "geometry/all.h"
 #include "partially_saturated_cells_method/all.h"

src/pe_coupling/discrete_particle_methods/all.h

+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla 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 3 of 
+//  the License, or (at your option) any later version.
+//  
+//  waLBerla 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 waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file all.h
+//! \ingroup pe_coupling
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+
+#include "waLBerlaDefinitions.h"
+
+#include "correlations/all.h"
+#include "evaluators/all.h"
+#include "gns_lbm/all.h"
+#include "utility/all.h"

src/pe_coupling/discrete_particle_methods/correlations/AddedMassForceCorrelations.h

+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla 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 3 of 
+//  the License, or (at your option) any later version.
+//  
+//  waLBerla 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 waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file AddedMassForceCorrelations.h
+//! \ingroup pe_coupling
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/DataTypes.h"
+
+#include <cmath>
+
+namespace walberla {
+namespace pe_coupling {
+namespace discrete_particle_methods {
+
+/*!\brief Correlation functions for the added mass force.
+ *
+ * To be compatible with the interface of AddedMassForceEvaluator, all functions use the following signature:
+ * Vector3<real_t> ( const Vector3<real_t> & timeDerivativeFluidVel, const Vector3<real_t> & timeDerivativeBodyVel,
+ *                   const real_t & bodyVolume, const real_t & fluidDensity )
+ */
+
+Vector3<real_t> addedMassForceFinn( const Vector3<real_t> & timeDerivativeFluidVel, const Vector3<real_t> & timeDerivativeBodyVel,
+                                    const real_t & bodyVolume, const real_t & fluidDensity )
+{
+   // formula from Finn et al(2016)
+   const real_t Coeffam = real_t(0.5);
+   return bodyVolume * fluidDensity * Coeffam * ( timeDerivativeFluidVel - timeDerivativeBodyVel );
+}
+
+Vector3<real_t> noAddedMassForce( const Vector3<real_t> &, const Vector3<real_t> &, const real_t &, const real_t & )
+{
+   return Vector3<real_t>(real_t(0));
+}
+
+} // namespace discrete_particle_methods
+} // namespace pe_coupling
+} // namespace walberla

src/pe_coupling/discrete_particle_methods/correlations/DragForceCorrelations.h

+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla 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 3 of 
+//  the License, or (at your option) any later version.
+//  
+//  waLBerla 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 waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file DragForceCorrelations.h
+//! \ingroup pe_coupling
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/math/Constants.h"
+
+namespace walberla {
+namespace pe_coupling {
+namespace discrete_particle_methods {
+
+
+/*!\brief Various correlation functions for the drag force.
+ *
+ * These functions calculate the drag force for fluid-particle interactions based on different empirical correlations
+ * from literature.
+ * Always be aware that those empirical formulas were obtained by different setups and are thus not generally applicable!
+ *
+ * To be compatible with the interface of InteractionForceEvaluator, all functions use the following signature:
+ * Vector3<real_t> ( const Vector3<real_t> & fluidVel, const Vector3<real_t> & particleVel, real_t solidVolumeFraction,
+ *                   real_t diameter, real_t fluidDynamicViscosity, real_t fluidDensity )
+ *
+ */
+
+
+// helper functions -> often used in more advanced drag correlations
+
+// equation to calculate the drag coefficient on isolated spherical particle
+// Schiller, L., Naumann, A., 1935. A drag coefficient correlation. Vdi Zeitung 77, 318-320.
+real_t dragCoeffSchillerNaumann( real_t reynoldsNumber )
+{
+   WALBERLA_ASSERT_GREATER_EQUAL( reynoldsNumber, real_t(0) );
+
+   return ( reynoldsNumber < real_t(1000) ) ? real_t(24) * ( real_t(1) + real_t(0.15) * std::pow(reynoldsNumber, real_t(0.687) ) ) / reynoldsNumber
+                                            : real_t(0.44);
+}
+
+// Coefficient from Stokes' law for drag, only valid for Stokes regime (low Reynolds numbers)
+// = 3 * PI * mu * D * fluidVolumeFraction
+real_t dragCoeffStokes ( real_t fluidVolumeFraction, real_t diameter, real_t fluidDynamicViscosity )
+{
+   return real_t(3) * math::M_PI * diameter * fluidDynamicViscosity * fluidVolumeFraction;
+}
+
+// threshold value for absolute relative velocity
+// if it is below this value, a drag force of 0 is set, to avoid instabilities stemming from divisions by this small value
+const real_t thresholdAbsoluteVelocityDifference = real_t(1e-10);
+
+
+//////////////////////
+//                  //
+//   CORRELATIONS   //
+//                  //
+//////////////////////
+
+// Stokes drag law
+Vector3<real_t> dragForceStokes( const Vector3<real_t> & fluidVel, const Vector3<real_t> & particleVel,
+                                 real_t solidVolumeFraction, real_t diameter, real_t fluidDynamicViscosity, real_t /*fluidDensity*/ )
+{
+   WALBERLA_ASSERT_GREATER_EQUAL( solidVolumeFraction, real_t(0) );
+   WALBERLA_ASSERT_LESS_EQUAL( solidVolumeFraction, real_t(1) );
+
+   Vector3<real_t> velDiff = fluidVel - particleVel;
+   real_t absVelDiff = velDiff.length();
+
+   if( absVelDiff < thresholdAbsoluteVelocityDifference ) return Vector3<real_t>(real_t(0));
+
+   real_t fluidVolumeFraction = real_t(1) - solidVolumeFraction;
+
+   return dragCoeffStokes( fluidVolumeFraction, diameter, fluidDynamicViscosity ) * velDiff;
+}
+
+
+// S. Ergun, Fluid flow through packed columns. Chemical Engineering Progress 48 (1952), 89-94.
+// Y. C. Wen, Y.H. Yu, Mechanics of fluidization. Chemical Engineering Progress Symposium Series 62 (1966), 100-111.
+// see also Beetstra, van der Hoef, Kuipers, "Drag Force of Intermediate Reynolds Number Flow Past Mono- and Bidisperse Arrays of Spheres" (2007)
+Vector3<real_t>  dragForceErgunWenYu( const Vector3<real_t> & fluidVel, const Vector3<real_t> & particleVel,
+                                      real_t solidVolumeFraction, real_t diameter, real_t fluidDynamicViscosity, real_t fluidDensity )
+{
+   WALBERLA_ASSERT_GREATER_EQUAL( solidVolumeFraction, real_t(0) );
+   WALBERLA_ASSERT_LESS_EQUAL( solidVolumeFraction, real_t(1) );
+
+   Vector3<real_t> velDiff = fluidVel - particleVel;
+   real_t absVelDiff = velDiff.length();
+
+   if( absVelDiff < thresholdAbsoluteVelocityDifference ) return Vector3<real_t>(real_t(0));
+
+   real_t fluidVolumeFraction = real_t(1) - solidVolumeFraction;
+
+   if( fluidVolumeFraction < real_t(0.8) )
+   {
+      // Ergun relation
+      real_t reynoldsNumber = fluidVolumeFraction * fluidDensity * absVelDiff * diameter / fluidDynamicViscosity;
+      real_t fDrag = real_t(150) * solidVolumeFraction / ( real_t(18) * fluidVolumeFraction * fluidVolumeFraction ) +
+                     real_t(1.75) / ( real_t(18) * fluidVolumeFraction * fluidVolumeFraction ) * reynoldsNumber;
+      return fDrag * dragCoeffStokes( fluidVolumeFraction, diameter, fluidDynamicViscosity ) * velDiff;
+   } else
+   {
+      // Wen & Yu correlation
+      real_t reynoldsNumber = fluidVolumeFraction * fluidDensity * absVelDiff * diameter / fluidDynamicViscosity;
+      real_t fDrag = dragCoeffSchillerNaumann( reynoldsNumber ) * reynoldsNumber / real_t(24) * std::pow( fluidVolumeFraction, real_t(-3.7) );
+      return fDrag * dragCoeffStokes( fluidVolumeFraction, diameter, fluidDynamicViscosity ) * velDiff;
+   }
+}
+
+// drag correlation proposed by Tang et al. - "A New Drag Correlation from Fully Resolved Simulations of Flow Past
+// Monodisperse Static Arrays of Spheres", AiChE, 2014
+Vector3<real_t> dragForceTang( const Vector3<real_t> & fluidVel, const Vector3<real_t> & particleVel,
+                               real_t solidVolumeFraction, real_t diameter, real_t fluidDynamicViscosity, real_t fluidDensity )
+{
+   WALBERLA_ASSERT_GREATER_EQUAL( solidVolumeFraction, real_t(0) );
+   WALBERLA_ASSERT_LESS_EQUAL( solidVolumeFraction, real_t(1) );
+
+   Vector3<real_t> velDiff = fluidVel - particleVel;
+   real_t absVelDiff = velDiff.length();
+
+   if( absVelDiff < thresholdAbsoluteVelocityDifference ) return Vector3<real_t>(real_t(0));
+
+   real_t fluidVolumeFraction = real_t(1) - solidVolumeFraction;
+   real_t fluidVolumeFractionP2 = fluidVolumeFraction * fluidVolumeFraction;
+   real_t inv_fluidVolumeFractionP4 = real_t(1) / (fluidVolumeFractionP2 * fluidVolumeFractionP2);
+   real_t reynoldsNumber = fluidVolumeFraction * fluidDensity * absVelDiff * diameter / fluidDynamicViscosity;
+
+   // Eq.21 from the paper
+   real_t fDrag = real_t(10) * solidVolumeFraction / fluidVolumeFractionP2 + fluidVolumeFractionP2 * ( real_t(1) + real_t(1.5) * std::sqrt(solidVolumeFraction) )
+                + ( real_t(0.11) * solidVolumeFraction * ( real_t(1) + solidVolumeFraction ) - real_t(0.00456) * inv_fluidVolumeFractionP4
+                + ( real_t(0.169) * fluidVolumeFraction + real_t(0.0644) * inv_fluidVolumeFractionP4 ) * std::pow( reynoldsNumber, -real_t(0.343) ) ) * reynoldsNumber;
+
+   return fDrag * dragCoeffStokes( fluidVolumeFraction, diameter, fluidDynamicViscosity ) * velDiff;
+
+}
+
+
+// drag correlation based on findings from Felice (1994)
+// used e.g. in Kafui et al (2002)
+Vector3<real_t> dragForceFelice( const Vector3<real_t> & fluidVel, const Vector3<real_t> & particleVel,
+                                 real_t solidVolumeFraction, real_t diameter, real_t fluidDynamicViscosity, real_t fluidDensity )
+{
+   WALBERLA_ASSERT_GREATER_EQUAL( solidVolumeFraction, real_t(0) );
+   WALBERLA_ASSERT_LESS_EQUAL( solidVolumeFraction, real_t(1) );
+
+   Vector3<real_t> velDiff = fluidVel - particleVel;
+   real_t absVelDiff = velDiff.length();
+
+   if( absVelDiff < thresholdAbsoluteVelocityDifference ) return Vector3<real_t>(real_t(0));
+
+   real_t fluidVolumeFraction = real_t(1) - solidVolumeFraction;
+
+   real_t reynoldsNumber = fluidVolumeFraction * fluidDensity * absVelDiff * diameter / fluidDynamicViscosity;
+
+   real_t temp1 = ( real_t(0.63) + real_t(4.8) / std::sqrt( reynoldsNumber ) );
+   real_t dragCoeff = temp1 * temp1;
+
+   real_t temp2 = real_t(1.5) - std::log10( reynoldsNumber );
+   real_t chi = real_t(3.7) - std::pow( real_t(0.65), (- real_t(0.5) * temp2 * temp2 ) );
+
+   return real_t(0.125) * dragCoeff * fluidDensity * math::M_PI * diameter * diameter * absVelDiff *
+          std::pow( fluidVolumeFraction, real_t(2) - chi) * velDiff;
+
+}
+
+// drag correlation based on findings from Tenneti, Garg, Subramaniam (2011)
+// used e.g. in Finn, Li, Apte - Particle based modelling and simulation of natural sand dynamics in the wave bottom boundary layer (2016)
+// could be generalized also for non-spherical particles, see Finn et al (2016)
+Vector3<real_t> dragForceTenneti( const Vector3<real_t> & fluidVel, const Vector3<real_t> & particleVel,
+                                  real_t solidVolumeFraction, real_t diameter, real_t fluidDynamicViscosity, real_t fluidDensity )
+{
+   WALBERLA_ASSERT_GREATER_EQUAL( solidVolumeFraction, real_t(0) );
+   WALBERLA_ASSERT_LESS_EQUAL( solidVolumeFraction, real_t(1) );
+
+   Vector3<real_t> velDiff = fluidVel - particleVel;
+   const real_t absVelDiff = velDiff.length();
+
+   if( absVelDiff < thresholdAbsoluteVelocityDifference ) return Vector3<real_t>(real_t(0));
+
+   const real_t fluidVolumeFraction = real_t(1) - solidVolumeFraction;
+
+   const real_t reynoldsNumber = fluidVolumeFraction * fluidDensity * absVelDiff * diameter / fluidDynamicViscosity;
+
+   const real_t fvfCubed = fluidVolumeFraction * fluidVolumeFraction * fluidVolumeFraction;
+   const real_t A = real_t(5.81) * solidVolumeFraction / fvfCubed + real_t(0.48) * std::cbrt( solidVolumeFraction ) / ( fvfCubed * fluidVolumeFraction );
+
+   const real_t svfCubed = solidVolumeFraction * solidVolumeFraction * solidVolumeFraction;
+   const real_t B = svfCubed * reynoldsNumber * ( real_t(0.95) + real_t(0.61) * svfCubed / ( fluidVolumeFraction * fluidVolumeFraction ) );
+
+   // version from Finn et al.
+   const real_t CdRe0Sphere = real_t(1) + real_t(0.15) *  std::pow( reynoldsNumber, real_t(0.687) );
+
+   const real_t CdRe = fluidVolumeFraction * ( CdRe0Sphere / fvfCubed + A + B );
+
+   return real_t(3) * math::M_PI * diameter * fluidDynamicViscosity * fluidVolumeFraction * CdRe * velDiff;
+
+}
+
+
+Vector3<real_t> noDragForce( const Vector3<real_t> & /*fluidVel*/, const Vector3<real_t> & /*particleVel*/,
+                             real_t /*solidVolumeFraction*/, real_t /*diameter*/, real_t /*fluidDynamicViscosity*/, real_t /*fluidDensity*/ )
+{
+   return Vector3<real_t>(real_t(0));
+}
+
+} // namespace discrete_particle_methods
+} // namespace pe_coupling
+} // namespace walberla

src/pe_coupling/discrete_particle_methods/correlations/LiftForceCorrelations.h

+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla 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 3 of 
+//  the License, or (at your option) any later version.
+//  
+//  waLBerla 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 waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file LiftForceCorrelations.h
+//! \ingroup pe_coupling
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/DataTypes.h"
+
+#include <cmath>
+
+namespace walberla {
+namespace pe_coupling {
+namespace discrete_particle_methods {
+
+/*!\brief Correlation functions for the lift force.
+ *
+ * To be compatible with the interface of LiftForceEvaluator, all functions use the following signature:
+ * Vector3<real_t> ( const Vector3<real_t> & fluidVel, const Vector3<real_t> & curlFluidVel, const Vector3<real_t> & particleVel,
+ *                   real_t diameter, real_t fluidDynamicViscosity, real_t fluidDensity )
+ */
+
+// Saffman lift force
+Vector3<real_t> liftForceSaffman ( const Vector3<real_t> & fluidVel, const Vector3<real_t> & curlFluidVel, const Vector3<real_t> & particleVel,
+                                   real_t diameter, real_t fluidDynamicViscosity, real_t fluidDensity )
+{
+   const real_t absCurlVel = curlFluidVel.length();
+   if( absCurlVel < real_t(1e-10) ) return Vector3<real_t>(real_t(0));
+
+   // Finn et al (2016) for spheres
+   const real_t Cl = real_t(1.61) * std::sqrt( ( fluidDynamicViscosity * fluidDensity) / absCurlVel );
+   return Cl * diameter * diameter * ( ( fluidVel - particleVel ) % curlFluidVel );
+
+   // Sun, Xiao (2016)
+   //const real_t Cl = real_t(1.6);
+   //return Cl * fluidDensity * std::sqrt( fluidDynamicViscosity / fluidDensity ) * diameter * diameter * ( ( fluidVel - particleVel ) % curlFluidVel );
+
+}
+
+Vector3<real_t> noLiftForce ( const Vector3<real_t> &, const Vector3<real_t> &, const Vector3<real_t> &, real_t, real_t, real_t )
+{
+   return Vector3<real_t>(real_t(0));
+}
+
+} // namespace discrete_particle_methods
+} // namespace pe_coupling
+} // namespace walberla

src/pe_coupling/discrete_particle_methods/correlations/all.h

+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla 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 3 of 
+//  the License, or (at your option) any later version.
+//  
+//  waLBerla 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 waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file all.h
+//! \ingroup pe_coupling
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+
+#include "waLBerlaDefinitions.h"
+
+#include "AddedMassForceCorrelations.h"
+#include "DragForceCorrelations.h"
+#include "LiftForceCorrelations.h"

src/pe_coupling/discrete_particle_methods/evaluators/AddedMassForceEvaluator.h

+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla 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 3 of 
+//  the License, or (at your option) any later version.
+//  
+//  waLBerla 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 waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file AddedMassForceEvaluator.h
+//! \ingroup pe_coupling
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/DataTypes.h"
+#include "core/debug/Debug.h"
+
+#include "domain_decomposition/StructuredBlockStorage.h"
+
+#include "field/interpolators/all.h"
+#include "field/distributors/all.h"
+#include "field/GhostLayerField.h"
+
+#include "lbm/field/MacroscopicValueCalculation.h"
+#include "lbm/field/PdfField.h"
+#include "lbm/lattice_model/ForceModel.h"
+
+#include "pe/rigidbody/BodyIterators.h"
+#include "pe/rigidbody/RigidBody.h"
+#include "pe/Types.h"
+#include "pe/Materials.h"
+
+#include "BodyVelocityTimeDerivativeEvaluator.h"
+
+#include "stencil/Directions.h"
+
+#include <boost/function.hpp>
+
+namespace walberla {
+namespace pe_coupling {
+namespace discrete_particle_methods {
+
+/*!\brief Evaluator of the added (virtual) mass force based on the given added mass correlation.
+ *
+ * The added mass force is evaluated for each body (requires interpolation of several quantities from the fluid fields)
+ * and then applied onto them.
+ * The corresponding reaction force is added to the fluid via the given distributor.
+ *
+ * Note that the fluid velocity (and its derivatives) has to be the fluid-phase velocity (and not the volume-averaged velocity).
+ *
+ * For more infos on interpolators, see field interpolators in src/field/interpolators.
+ * For more infos on distributors, see src/field/distributors.
+ */
+template< typename FlagField_T, template<typename,typename> class FieldInterpolator_T, template<typename,typename> class Distributor_T >
+class AddedMassForceEvaluator
+{  
+
+public:
+
+   typedef GhostLayerField< Vector3<real_t>, 1>             Vec3Field_T;
+   typedef FieldInterpolator_T<Vec3Field_T, FlagField_T>    Vec3FieldInterpolator_T;
+   typedef Distributor_T<Vec3Field_T, FlagField_T>          ForceDistributor_T;
+
+   AddedMassForceEvaluator( const shared_ptr<StructuredBlockStorage> & blockStorage,
+                            const BlockDataID & forceFieldID, const BlockDataID & bodyStorageID,
+                            const BlockDataID & flagFieldID, const Set< FlagUID > & domainFlags,
+                            const BlockDataID & velocityTimeDerivativeFieldID,
+                            const boost::function<Vector3<real_t> ( const Vector3<real_t> &, const Vector3<real_t> &, real_t, real_t )> & addedMassForceCorrelationFunction,
+                            const shared_ptr< BodyVelocityTimeDerivativeEvaluator > & bodyVelocityTimeDerivativeEvaluator )
+   : blockStorage_( blockStorage ), bodyStorageID_( bodyStorageID ),
+     addedMassForceCorrelationFunction_( addedMassForceCorrelationFunction ), bodyVelocityTimeDerivativeEvaluator_( bodyVelocityTimeDerivativeEvaluator )
+   {
+      velocityTimeDerivativeFieldInterpolatorID_ = field::addFieldInterpolator< Vec3FieldInterpolator_T, FlagField_T >( blockStorage, velocityTimeDerivativeFieldID, flagFieldID, domainFlags );
+      forceDistributorID_ = field::addDistributor< ForceDistributor_T, FlagField_T >( blockStorage, forceFieldID, flagFieldID, domainFlags );
+   }
+
+   void operator()();
+
+private:
+
+   shared_ptr<StructuredBlockStorage> blockStorage_;
+   const BlockDataID bodyStorageID_;
+
+   boost::function<Vector3<real_t> ( const Vector3<real_t> &, const Vector3<real_t> &, real_t, real_t )> addedMassForceCorrelationFunction_;
+
+   shared_ptr< BodyVelocityTimeDerivativeEvaluator > bodyVelocityTimeDerivativeEvaluator_;
+
+   BlockDataID velocityTimeDerivativeFieldInterpolatorID_;
+   BlockDataID forceDistributorID_;
+
+};
+
+template< typename FlagField_T, template<typename,typename> class FieldInterpolator_T, template<typename,typename> class Distributor_T >
+void AddedMassForceEvaluator< FlagField_T, FieldInterpolator_T, Distributor_T >
+::operator()()
+{
+
+   for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
+   {
+
+      Vec3FieldInterpolator_T * velocityTimeDerivativeInterpolator = blockIt->getData< Vec3FieldInterpolator_T >( velocityTimeDerivativeFieldInterpolatorID_ );
+      ForceDistributor_T * forceDistributor                        = blockIt->getData< ForceDistributor_T >( forceDistributorID_ );
+
+      for( auto bodyIt = pe::LocalBodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt )
+      {
+         //TODO check if body should be treated by DPM
+
+         Vector3<real_t> forceOnFluid( real_t(0) );
+
+         Vector3<real_t> bodyPosition = bodyIt->getPosition();
+         Vector3<real_t> bodyVelocity = bodyIt->getLinearVel();
+         real_t bodyVolume = bodyIt->getVolume();
+         real_t fluidDensity( real_t(1) );
+
+         // evaluate added (virtual) mass force
+         Vector3<real_t> timeDerivativeFluidVelocity( real_t(0) );
+         Vector3<real_t> timeDerivativeBodyVelocity( real_t(0) );
+         bodyVelocityTimeDerivativeEvaluator_->get(timeDerivativeBodyVelocity, bodyVelocity, bodyIt->getSystemID() );
+         velocityTimeDerivativeInterpolator->get( bodyPosition, &timeDerivativeFluidVelocity );
+         Vector3<real_t> addedMassForce = addedMassForceCorrelationFunction_( timeDerivativeFluidVelocity, timeDerivativeBodyVelocity, bodyVolume, fluidDensity );
+         WALBERLA_ASSERT( !math::isnan(addedMassForce[0]) && !math::isnan(addedMassForce[1]) && !math::isnan(addedMassForce[2]),
+                          "NaN found in added mass force for body at position " << bodyPosition );
+
+         bodyIt->addForce( addedMassForce );
+         forceOnFluid += ( -addedMassForce );
+
+         // set/distribute force on fluid
+         forceDistributor->distribute( bodyPosition, &forceOnFluid );
+      }
+   }
+
+}
+
+} // namespace discrete_particle_methods
+} // namespace pe_coupling
+} // namespace walberla