path: root/examples/particles/bifurcation3d/eulerEuler/bifurcation3d.cpp

/*  Lattice Boltzmann sample, written in C++, using the OpenLB
 *  library
 *
 *  Copyright (C) 2011-2016 Robin Trunk, Mathias J. Krause
 *  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.
 */

/* bifurcation3d.cpp:
 * This example examines a steady particulate flow past a bifurcation. At the inlet,
 * an inflow condition with grad_n u = 0 and rho = 1 is implemented.
 * At both outlets, a Poiseuille profile is imposed on the velocity.
 * After a start time, particles are put into the bifurcation by imposing
 * a inflow condition with rho = 1 on the second euler phase at the inlet.
 * The particles are simulated as a continuum with a advection-diffusion equation
 * and experience a stokes drag force.
 *
 * A publication using the same geometry can be found here:
 * http://link.springer.com/chapter/10.1007/978-3-642-36961-2_5
 *  *
 */

#include "olb3D.h"
#include "olb3D.hh"   // use only generic version!

using namespace std;
using namespace olb;
using namespace olb::descriptors;
using namespace olb::graphics;
using namespace olb::util;

typedef double T;
#define NSDESCRIPTOR D3Q19<>
#define ADDESCRIPTOR D3Q7<VELOCITY,VELOCITY2>

const T Re = 50;               // Reynolds number
const int N = 19;              // resolution of the model
const int iTperiod = 100;      // amount of timesteps when new boundary conditions are reset and results are visualized
const T diffusion = 1.e-6;     // diffusion coefficient for advection-diffusion equation
const T radius = 1.5e-04;      // particles radius
const T partRho = 998.2;       // particles density
const T maxPhysT = 10.;        // max. simulation time in s, SI unit

// center of inflow and outflow regions [m]
Vector<T,3> inletCenter( T(), T(), 0.0786395 );
Vector<T,3> outletCenter0( -0.0235929682287551, -0.000052820468762797, -0.021445708949909 );
Vector<T,3> outletCenter1( 0.0233643529416147, 0.00000212439067050152, -0.0211994104877918 );

// radii of inflow and outflow regions [m]
T inletRadius = 0.00999839;
T outletRadius0 = 0.007927;
T outletRadius1 = 0.00787134;

// normals of inflow and outflow regions
Vector<T,3> inletNormal( T(), T(), T( -1 ) );
Vector<T,3> outletNormal0( 0.505126, -0.04177, 0.862034 );
Vector<T,3> outletNormal1( -0.483331, -0.0102764, 0.875377 );

void prepareGeometry( UnitConverter<T,NSDESCRIPTOR> const& converter,
                      IndicatorF3D<T>& indicator, STLreader<T>& stlReader,
                      SuperGeometry3D<T>& superGeometry )
{

  OstreamManager clout( std::cout, "prepareGeometry" );

  clout << "Prepare Geometry ..." << std::endl;

  superGeometry.rename( 0, 2, indicator );
  superGeometry.rename( 2, 1, stlReader );

  superGeometry.clean();

  // rename the material at the inlet
  IndicatorCircle3D<T> inletCircle( inletCenter, inletNormal, inletRadius );
  IndicatorCylinder3D<T> inlet( inletCircle, 2 * converter.getConversionFactorLength() );
  superGeometry.rename( 2, 3, 1, inlet );

  // rename the material at the outlet0
  IndicatorCircle3D<T> outletCircle0( outletCenter0, outletNormal0, 0.95*outletRadius0 );
  IndicatorCylinder3D<T> outlet0( outletCircle0, 4 * converter.getConversionFactorLength() );
  superGeometry.rename( 2, 4, outlet0 );

  // rename the material at the outlet1
  IndicatorCircle3D<T> outletCircle1( outletCenter1, outletNormal1, 0.95*outletRadius1 );
  IndicatorCylinder3D<T> outlet1( outletCircle1, 4 * converter.getConversionFactorLength() );
  superGeometry.rename( 2, 5, outlet1 );

  superGeometry.clean();
  superGeometry.innerClean( true );
  superGeometry.checkForErrors();

  superGeometry.print();

  clout << "Prepare Geometry ... OK" << std::endl;
  return;
}

void prepareLattice( SuperLattice3D<T, NSDESCRIPTOR>& sLatticeNS,
                     SuperLattice3D<T, ADDESCRIPTOR>& sLatticeAD,
                     UnitConverter<T,NSDESCRIPTOR> const& converter,
                     Dynamics<T, NSDESCRIPTOR>& bulkDynamics,
                     Dynamics<T, ADDESCRIPTOR>& bulkDynamicsAD,
                     sOnLatticeBoundaryCondition3D<T, NSDESCRIPTOR>& bc,
                     sOnLatticeBoundaryCondition3D<T, ADDESCRIPTOR>& bcAD,
                     SuperGeometry3D<T>& superGeometry,
                     T omegaAD )
{

  OstreamManager clout( std::cout, "prepareLattice" );
  clout << "Prepare Lattice ..." << std::endl;

  const T omega = converter.getLatticeRelaxationFrequency();

  // Material=0 --> do nothing
  sLatticeNS.defineDynamics( superGeometry, 0, &instances::getNoDynamics<T, NSDESCRIPTOR>() );
  sLatticeAD.defineDynamics( superGeometry, 0, &instances::getNoDynamics<T, ADDESCRIPTOR>() );

  // Material=1 --> bulk dynamics
  // Material=3 --> bulk dynamics (inflow)
  auto inflowIndicator = superGeometry.getMaterialIndicator({1, 3});
  sLatticeNS.defineDynamics( inflowIndicator, &bulkDynamics );
  sLatticeAD.defineDynamics( inflowIndicator, &bulkDynamicsAD );

  // Material=2 --> bounce-back /
  sLatticeNS.defineDynamics( superGeometry, 2, &instances::getBounceBack<T, NSDESCRIPTOR>() );
  sLatticeAD.defineDynamics( superGeometry, 2, &instances::getZeroDistributionDynamics<T, ADDESCRIPTOR>() );

  // Material=4,5 -->bulk dynamics / do-nothing (outflow)
  auto outflowIndicator = superGeometry.getMaterialIndicator({4, 5});
  sLatticeNS.defineDynamics( outflowIndicator, &bulkDynamics );
  sLatticeAD.defineDynamics( outflowIndicator, &instances::getNoDynamics<T, ADDESCRIPTOR>() );

  // Setting of the boundary conditions
  bc.addPressureBoundary( superGeometry, 3, omega );
  bc.addVelocityBoundary( outflowIndicator, omega );
  bcAD.addZeroDistributionBoundary( superGeometry, 2 );
  bcAD.addTemperatureBoundary( superGeometry, 3, omegaAD );
  bcAD.addConvectionBoundary( outflowIndicator );
  bcAD.addExtFieldBoundary( superGeometry.getMaterialIndicator({2, 3, 4, 5}), ADDESCRIPTOR::index<descriptors::VELOCITY>() );

  // Initial conditions
  AnalyticalConst3D<T,T> rho1( 1. );
  AnalyticalConst3D<T,T> rho0( 1.e-8 );
  std::vector<T> velocity( 3,T() );
  AnalyticalConst3D<T,T> u0( velocity );

  // Initialize all values of distribution functions to their local equilibrium
  sLatticeNS.defineRhoU( superGeometry.getMaterialIndicator({1, 2, 3, 4, 5}), rho1, u0 );
  sLatticeNS.iniEquilibrium( superGeometry.getMaterialIndicator({1, 2, 3, 4, 5}), rho1, u0 );
  sLatticeAD.defineRho( superGeometry, 3, rho1 );
  sLatticeAD.iniEquilibrium( superGeometry.getMaterialIndicator({1, 2, 4, 5}), rho0, u0 );

  // Lattice initialize
  sLatticeNS.initialize();
  sLatticeAD.initialize();

  clout << "Prepare Lattice ... OK" << std::endl;
  return;
}

void setBoundaryValues( SuperLattice3D<T, NSDESCRIPTOR>& sLatticeNS,
                        UnitConverter<T,NSDESCRIPTOR> const& converter, int iT,
                        T maxPhysT, SuperGeometry3D<T>& superGeometry )
{

  OstreamManager clout( std::cout, "setBoundaryValues" );

  // No of time steps for smooth start-up
  int iTmaxStart = converter.getLatticeTime( 0.8*maxPhysT );
  // Set inflow velocity
  T maxVelocity = converter.getCharLatticeVelocity() * 3. / 4. * std::pow(
      inletRadius, 2 ) / std::pow( outletRadius0, 2 );
  if ( iT % iTperiod == 0 ) {
    if ( iT <= iTmaxStart ) {
      SinusStartScale<T, int> startScale( iTmaxStart,