flow_past_cylinder.hpp

Go to the documentation of this file.

#pragma once
 
#include <deal.II/base/function.h>
#include <deal.II/base/function_signed_distance.h>
#include <deal.II/base/point.h>
#include <deal.II/base/tensor_function.h>
 
#include <deal.II/distributed/tria.h>
 
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/manifold_lib.h>
 
#include <meltpooldg/core/simulation_case_base.hpp>
#include <meltpooldg/utilities/characteristic_functions.hpp>
#include <meltpooldg/utilities/utility_functions.hpp>
 
#include <iostream>
 
#include "../../melt_pool_case.hpp"
 
namespace MeltPoolDG::Simulation::FlowPastCylinder
{
  template <int dim>
  class InitializePhi : public dealii::Function<dim>
  {
  public:
    InitializePhi()
      : dealii::Function<dim>()
      , distance_sphere(dim == 1 ? dealii::Point<dim>(0.0) : dealii::Point<dim>(1.3, 0.2), 0.1)
    {}
 
    double
    value(const dealii::Point<dim> &p, const unsigned int /*component*/) const override
    {
      return CharacteristicFunctions::sgn(-distance_sphere.value(p));
    }
 
  private:
    const dealii::Functions::SignedDistance::Sphere<dim> distance_sphere;
  };
 
  template <int dim>
  class InflowVelocity : public dealii::Function<dim>
  {
  public:
    InflowVelocity(const double time, const bool fluctuating)
      : dealii::Function<dim>(dim, time)
      , fluctuating(fluctuating)
    {}
 
    void
    vector_value(const dealii::Point<dim> &p, dealii::Vector<double> &values) const override
    {
      AssertDimension(values.size(), dim);
 
      // inflow velocity according to Schaefer & Turek
      const double Um          = (dim == 2 ? 1.5 : 2.25);
      const double H           = 0.41;
      double       coefficient = dealii::Utilities::fixed_power<dim - 1>(4.) * Um /
                           dealii::Utilities::fixed_power<2 * dim - 2>(H);
      values(0) = coefficient * p[1] * (H - p[1]);
      if (dim == 3)
        values(0) *= p[2] * (H - p[2]);
      if (fluctuating)
        values(0) *= std::sin(this->get_time() * dealii::numbers::PI / 8.);
      for (unsigned int d = 1; d < dim; ++d)
        values(d) = 0;
    }
 
  private:
    const bool fluctuating;
  };
 
  /* for constant Dirichlet conditions we could also use the ConstantFunction
   * utility from dealii
   */
  template <int dim>
  class DirichletCondition : public dealii::Function<dim>
  {
  public:
    DirichletCondition()
      : dealii::Function<dim>()
    {}
 
    double
    value(const dealii::Point<dim> &p, const unsigned int component = 0) const override
    {
      (void)p;
      (void)component;
      return -1.0;
    }
  };
 
  /*
   *      This class collects all relevant input data for the level set simulation
   */
 
  template <int dim, typename number>
  class SimulationFlowPastCylinder : public MeltPoolCase<dim, number>
  {
  public:
    SimulationFlowPastCylinder(std::string parameter_file, const MPI_Comm mpi_communicator)
      : MeltPoolCase<dim, number>(parameter_file, mpi_communicator)
    {}
 
    void
    create_spatial_discretization() override
    {
      this->triangulation =
        std::make_shared<dealii::parallel::distributed::Triangulation<dim>>(this->mpi_communicator);
 
      if constexpr (dim == 2)
        {
          dealii::Point<2>          p1(0, 0);
          dealii::Point<2>          p2(2.5, 0.4);
          std::vector<unsigned int> refinements({50, 8});
          dealii::Triangulation<2>  tmp;
          dealii::GridGenerator::subdivided_hyper_rectangle(tmp, refinements, p1, p2);
          std::set<dealii::Triangulation<2>::active_cell_iterator> cells_in_void;
          for (dealii::Triangulation<2>::active_cell_iterator cell = tmp.begin(); cell != tmp.end();
               ++cell)
            if (cell->center()[0] > 0.45 && cell->center()[0] < 0.55 && cell->center()[1] > 0.15 &&
                cell->center()[1] < 0.25)
              cells_in_void.insert(cell);
          dealii::GridGenerator::create_triangulation_with_removed_cells(tmp,
                                                                         cells_in_void,
                                                                         *this->triangulation);
 
          // shift cells at the upper end of the domain from 0.40 to 0.41. It
          // corresponds to faces with id 3
          for (dealii::Triangulation<2>::cell_iterator cell = this->triangulation->begin();
               cell != this->triangulation->end();
               ++cell)
            if (cell->at_boundary(3) && cell->face(3)->center()[1] > 0.39999999999)
              for (unsigned int v = 0; v < dealii::GeometryInfo<2>::vertices_per_face; ++v)
                cell->face(3)->vertex(v)[1] = 0.41;
 
          // Set the left boundary (inflow) to 1, the right to 2, the rest to 0.
          for (dealii::Triangulation<2>::active_cell_iterator cell = this->triangulation->begin();
               cell != this->triangulation->end();
               ++cell)
            for (unsigned int f = 0; f < dealii::GeometryInfo<2>::faces_per_cell; ++f)
              if (cell->face(f)->at_boundary())
                {
                  if (std::abs(cell->face(f)->center()[0]) < 1e-12)
                    cell->face(f)->set_all_boundary_ids(1);
                  else if (std::abs(cell->face(f)->center()[0] - 2.5) < 1e-12)
                    cell->face(f)->set_all_boundary_ids(2);
                  else
                    cell->face(f)->set_all_boundary_ids(0);
                }
        }
      else
        {
          AssertThrow(false, dealii::ExcNotImplemented());
        }
    }
 
    void
    set_boundary_conditions() override
    {
      auto dirichlet = std::make_shared<DirichletCondition<dim>>();
 
      this->attach_boundary_condition({0, dirichlet}, "dirichlet", "level_set");
      this->attach_boundary_condition({1, dirichlet}, "dirichlet", "level_set");
 
      this->attach_boundary_condition(0, "no_slip", "navier_stokes_u");
 
      this->attach_boundary_condition(
        {1, std::shared_ptr<dealii::Function<dim>>(new InflowVelocity<dim>(0., true))},
        "dirichlet",
        "navier_stokes_u");
 
      this->attach_boundary_condition(
        {2, std::shared_ptr<dealii::Function<dim>>(new dealii::Functions::ZeroFunction<dim>(1))},
        "neumann",
        "navier_stokes_p");
    }
 
    void
    set_field_conditions() override
    {
      this->attach_initial_condition(std::make_shared<InitializePhi<dim>>(), "level_set");
      this->attach_initial_condition(std::shared_ptr<dealii::Function<dim>>(
                                       new InflowVelocity<dim>(0., false)),
                                     "navier_stokes_u");
    }
  };
} // namespace MeltPoolDG::Simulation::FlowPastCylinder