Tracer drift flux model (#1014)

Description of the problem
For some simple cases of tracer modeling, we need to account for tracer drift to have a simple model for dilute multiphase simulation. This is currently not possible.

Description of the solution
Add a simple drift flux velocity model for the tracer physics. This is just a user-defined function that adds a supplementary velocity on top of the fluid velocity to the tracer physics.

How Has This Been Tested?
This is a very minimal modification which just adds a velocity to the fluid velocity so this does not require specific unit test. If all the regular tests still pass, this implementation should be valid.

Documentation
Documentation of the class was added and class is adequately documented.

doc/source/parameters/cfd/cfd.rst

@@ -31,5 +31,6 @@ General, CFD and Multiphysics
    source_term
    stabilization
    timer
-   velocity_source.rst
+   tracer_drift_velocity
+   velocity_source
    volume_of_fluid

doc/source/parameters/cfd/tracer_drift_velocity.rst

@@ -0,0 +1,15 @@
+=====================
+Tracer Drift Velocity
+=====================
+
+This subsection allows you to define a tracer drift velocity. This drift velocity is an additional velocity which is added to the fluid velocity when advecting the tracer. This enables user to model, with some strong hypothesis, the dynamics of a dispersed phase (such as bubbles or particles) within a fluid flow. 
+
+.. code-block:: text
+
+  subsection tracer drift velocity
+    subsection drift velocity
+      # Default values in 2D
+      set Function expression = 0; 0
+      # in 3D: set Function expression = 0; 0; 0
+    end
+  end

include/solvers/analytical_solutions.h

@@ -87,7 +87,7 @@ namespace AnalyticalSolutions
     declare_parameters(ParameterHandler &prm);
 
     /**
-     * @brief Declares the parameters required by the analytical solution
+     * @brief Parses the parameters required by the analytical solution
      * within a parameter file
      *
      * @param prm ParameterHandler used to parse the parameters.

include/solvers/simulation_parameters.h

@@ -32,6 +32,7 @@
 #include <solvers/initial_conditions.h>
 #include <solvers/physical_properties_manager.h>
 #include <solvers/source_terms.h>
+#include <solvers/tracer_drift_velocity.h>
 
 template <int dim>
 class SimulationParameters
@@ -69,6 +70,7 @@ class SimulationParameters
   Parameters::Stabilization                     stabilization;
   Parameters::ALE<dim>                          ale;
   Parameters::Evaporation                       evaporation;
+  Parameters::TracerDriftVelocity<dim>          tracer_drift_velocity;
 
 
 
@@ -138,6 +140,8 @@ class SimulationParameters
     Parameters::Evaporation::declare_parameters(prm);
 
     multiphysics.declare_parameters(prm);
+
+    tracer_drift_velocity.declare_parameters(prm);
   }
 
   void
@@ -182,6 +186,7 @@ class SimulationParameters
     stabilization.parse_parameters(prm);
     ale.parse_parameters(prm);
     evaporation.parse_parameters(prm);
+    tracer_drift_velocity.parse_parameters(prm);
 
     physical_properties_manager.initialize(physical_properties);
 

include/solvers/tracer_drift_velocity.h

@@ -0,0 +1,108 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2019 -  by the Lethe authors
+ *
+ * This file is part of the Lethe library
+ *
+ * The Lethe library is free software; you can use it, redistribute
+ * it, and/or modify it under the terms of the GNU Lesser General
+ * Public License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ * The full text of the license can be found in the file LICENSE at
+ * the top level of the Lethe distribution.
+ *
+ * ---------------------------------------------------------------------
+ *
+ */
+
+
+#ifndef lethe_tracer_drift_velocity_h
+#define lethe_tracer_drift_velocity_h
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/parameter_handler.h>
+#include <deal.II/base/parsed_function.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <memory>
+
+using namespace dealii;
+/**
+ * Drift velocity that is added to the tracer velocity to account for simplified
+ *multiphase simulations using drift-flux modeling.
+ **/
+
+namespace Parameters
+{
+  /**
+* @brief Implements a drift velocity function to account for simplified multiphase flows.
+*
+* @tparam dim An integer that denotes the dimension of the space in which
+* the flow is solved.
+
+  * The drift velocity provides a simple way to model dilute disperse multiphase
+flow through the tracer physics.
+  **/
+
+  template <int dim>
+  class TracerDriftVelocity
+  {
+  public:
+    TracerDriftVelocity()
+    {
+      drift_velocity = std::make_shared<Functions::ParsedFunction<dim>>(dim);
+    }
+
+    /**
+     * @brief Declares the parameters required by the drift velocity.
+     *
+     * @param prm ParameterHandler used to declare the parameters.
+     *
+     */
+    virtual void
+    declare_parameters(ParameterHandler &prm);
+
+    /**
+     * @brief Parses the parameters required by the analytical solution
+     * within a parameter file
+     *
+     * @param prm ParameterHandler used to parse the parameters.
+     *
+     */
+    virtual void
+    parse_parameters(ParameterHandler &prm);
+
+    // Drift velocity
+    std::shared_ptr<Functions::ParsedFunction<dim>> drift_velocity;
+  };
+
+  template <int dim>
+  void
+  TracerDriftVelocity<dim>::declare_parameters(ParameterHandler &prm)
+  {
+    prm.enter_subsection("tracer drift velocity");
+
+    prm.enter_subsection("drift velocity");
+    drift_velocity->declare_parameters(prm, dim);
+    prm.leave_subsection();
+
+    prm.leave_subsection();
+  }
+
+  template <int dim>
+  void
+  TracerDriftVelocity<dim>::parse_parameters(ParameterHandler &prm)
+  {
+    prm.enter_subsection("tracer drift velocity");
+
+    prm.enter_subsection("drift velocity");
+    drift_velocity->parse_parameters(prm);
+    prm.leave_subsection();
+
+    prm.leave_subsection();
+  }
+} // namespace Parameters
+
+
+#endif

include/solvers/tracer_scratch_data.h

@@ -208,19 +208,38 @@ class TracerScratchData
    *
    * @param ale The ALE parameters which include the ALE function
    *
+   * @param drift_velocity Function to calculate the drift velocity
+   *
    */
 
   template <typename VectorType>
   void
-  reinit_velocity(const typename DoFHandler<dim>::active_cell_iterator &cell,
-                  const VectorType           &current_solution,
-                  const Parameters::ALE<dim> &ale)
+  reinit_velocity(
+    const typename DoFHandler<dim>::active_cell_iterator &cell,
+    const VectorType                                     &current_solution,
+    const Parameters::ALE<dim>                           &ale,
+    std::shared_ptr<Functions::ParsedFunction<dim>>       drift_velocity)
   {
     this->fe_values_fd.reinit(cell);
 
     this->fe_values_fd[velocities].get_function_values(current_solution,
                                                        velocity_values);
 
+    // Add the drift velocity to the velocity to account for tracer drift flux
+    // modeling
+    Tensor<1, dim> drift_velocity_tensor;
+    Vector<double> drift_velocity_vector(dim);
+
+    for (unsigned int q = 0; q < n_q_points; ++q)
+      {
+        drift_velocity->vector_value(quadrature_points[q],
+                                     drift_velocity_vector);
+        for (unsigned int d = 0; d < dim; ++d)
+          drift_velocity_tensor[d] = drift_velocity_vector[d];
+
+        velocity_values[q] += drift_velocity_tensor;
+      }
+
     if (!ale.enabled())
       return;
 

source/solvers/CMakeLists.txt

@@ -31,6 +31,7 @@ add_library(lethe-solvers
   source_terms.cc
   tracer.cc
   tracer_assemblers.cc
+  tracer_drift_velocity.cc
   tracer_scratch_data.cc
   vof.cc
   vof_assemblers.cc
@@ -71,6 +72,7 @@ add_library(lethe-solvers
   ../../include/solvers/stabilization.h
   ../../include/solvers/tracer.h
   ../../include/solvers/tracer_assemblers.h
+  ../../include/solvers/tracer_drift_velocity.h
   ../../include/solvers/tracer_scratch_data.h
   ../../include/solvers/vof.h
   ../../include/solvers/vof_assemblers.h

source/solvers/tracer.cc

@@ -109,14 +109,16 @@ Tracer<dim>::assemble_local_system_matrix(
             velocity_cell,
             *multiphysics->get_block_time_average_solution(
               PhysicsID::fluid_dynamics),
-            this->simulation_parameters.ale);
+            this->simulation_parameters.ale,
+            this->simulation_parameters.tracer_drift_velocity.drift_velocity);
         }
       else
         {
-          scratch_data.reinit_velocity(velocity_cell,
-                                       *multiphysics->get_block_solution(
-                                         PhysicsID::fluid_dynamics),
-                                       this->simulation_parameters.ale);
+          scratch_data.reinit_velocity(
+            velocity_cell,
+            *multiphysics->get_block_solution(PhysicsID::fluid_dynamics),
+            this->simulation_parameters.ale,
+            this->simulation_parameters.tracer_drift_velocity.drift_velocity);
         }
     }
   else
@@ -126,17 +128,19 @@ Tracer<dim>::assemble_local_system_matrix(
           simulation_control->get_current_time() >
             this->simulation_parameters.post_processing.initial_time)
         {
-          scratch_data.reinit_velocity(velocity_cell,
-                                       *multiphysics->get_time_average_solution(
-                                         PhysicsID::fluid_dynamics),
-                                       this->simulation_parameters.ale);
+          scratch_data.reinit_velocity(
+            velocity_cell,
+            *multiphysics->get_time_average_solution(PhysicsID::fluid_dynamics),
+            this->simulation_parameters.ale,
+            this->simulation_parameters.tracer_drift_velocity.drift_velocity);
         }
       else
         {
-          scratch_data.reinit_velocity(velocity_cell,
-                                       *multiphysics->get_solution(
-                                         PhysicsID::fluid_dynamics),
-                                       this->simulation_parameters.ale);
+          scratch_data.reinit_velocity(
+            velocity_cell,
+            *multiphysics->get_solution(PhysicsID::fluid_dynamics),
+            this->simulation_parameters.ale,
+            this->simulation_parameters.tracer_drift_velocity.drift_velocity);
         }
     }
 
@@ -235,14 +239,16 @@ Tracer<dim>::assemble_local_system_rhs(
             velocity_cell,
             *multiphysics->get_block_time_average_solution(
               PhysicsID::fluid_dynamics),
-            this->simulation_parameters.ale);
+            this->simulation_parameters.ale,
+            this->simulation_parameters.tracer_drift_velocity.drift_velocity);
         }
       else
         {
-          scratch_data.reinit_velocity(velocity_cell,
-                                       *multiphysics->get_block_solution(
-                                         PhysicsID::fluid_dynamics),
-                                       this->simulation_parameters.ale);
+          scratch_data.reinit_velocity(
+            velocity_cell,
+            *multiphysics->get_block_solution(PhysicsID::fluid_dynamics),
+            this->simulation_parameters.ale,
+            this->simulation_parameters.tracer_drift_velocity.drift_velocity);
         }
     }
   else
@@ -252,17 +258,19 @@ Tracer<dim>::assemble_local_system_rhs(
           simulation_control->get_current_time() >
             this->simulation_parameters.post_processing.initial_time)
         {
-          scratch_data.reinit_velocity(velocity_cell,
-                                       *multiphysics->get_time_average_solution(
-                                         PhysicsID::fluid_dynamics),
-                                       this->simulation_parameters.ale);
+          scratch_data.reinit_velocity(
+            velocity_cell,
+            *multiphysics->get_time_average_solution(PhysicsID::fluid_dynamics),
+            this->simulation_parameters.ale,
+            this->simulation_parameters.tracer_drift_velocity.drift_velocity);
         }
       else
         {
-          scratch_data.reinit_velocity(velocity_cell,
-                                       *multiphysics->get_solution(
-                                         PhysicsID::fluid_dynamics),
-                                       this->simulation_parameters.ale);
+          scratch_data.reinit_velocity(
+            velocity_cell,
+            *multiphysics->get_solution(PhysicsID::fluid_dynamics),
+            this->simulation_parameters.ale,
+            this->simulation_parameters.tracer_drift_velocity.drift_velocity);
         }
     }
 

source/solvers/tracer_drift_velocity.cc

@@ -0,0 +1,22 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2019 - by the Lethe authors
+ *
+ * This file is part of the Lethe library
+ *
+ * The Lethe library is free software; you can use it, redistribute
+ * it, and/or modify it under the terms of the GNU Lesser General
+ * Public License as published by the Free Software Foundation; either
+ * version 3.1 of the License, or (at your option) any later version.
+ * The full text of the license can be found in the file LICENSE at
+ * the top level of the Lethe distribution.
+ *
+ * ---------------------------------------------------------------------
+ */
+
+#include <solvers/tracer_drift_velocity.h>
+
+// Pre-compile the 2D and 3D Navier-Stokes solver to ensure that the library is
+// valid before we actually compile the solver This greatly helps with debugging
+template class Parameters::TracerDriftVelocity<2>;
+template class Parameters::TracerDriftVelocity<3>;