hc process revision

Documentation/ProjectFile/material/fluid/density/ConcentrationAndPressureDependent/c_ConcentrationAndPressureDependent.md

@@ -0,0 +1 @@
+\copydoc MaterialLib::Fluid::LinearConcentrationAndPressureDependentDensity

Documentation/ProjectFile/material/fluid/density/ConcentrationAndPressureDependent/t_fluid_density_concentration_difference_ratio.md

@@ -0,0 +1 @@
+\f$ \bar \alpha\f$ is difference ratio for fluid density depending on concentration

Documentation/ProjectFile/material/fluid/density/ConcentrationAndPressureDependent/t_fluid_density_pressure_difference_ratio.md

@@ -0,0 +1 @@
+\f$ \bar \alpha\f$ is difference ratio for fluid density depending on pressure

Documentation/ProjectFile/material/fluid/density/ConcentrationAndPressureDependent/t_reference_concentration.md

@@ -0,0 +1 @@
+\f$ C_{\text{ref}}\f$ is the reference concentration

Documentation/ProjectFile/material/fluid/density/ConcentrationAndPressureDependent/t_reference_density.md

@@ -0,0 +1 @@
+\f$ \varrho_{\text{ref}}\f$ is the reference density

Documentation/ProjectFile/material/fluid/density/ConcentrationAndPressureDependent/t_reference_pressure.md

@@ -0,0 +1 @@
+\f$ p_{\text{ref}}\f$ is the reference pressure

MaterialLib/Fluid/Density/CreateFluidDensityModel.cpp

@@ -16,6 +16,7 @@
 #include "BaseLib/Error.h"
 
 #include "IdealGasLaw.h"
+#include "LinearConcentrationAndPressureDependentDensity.h"
 #include "LinearConcentrationDependentDensity.h"
 #include "LinearTemperatureDependentDensity.h"
 #include "LiquidDensity.h"
@@ -91,6 +92,37 @@ static std::unique_ptr<FluidProperty> createLinearConcentrationDependentDensity(
         reference_concentration,
         fluid_density_difference_ratio);
 }
+static std::unique_ptr<FluidProperty>
+createLinearConcentrationAndPressureDependentDensity(
+    BaseLib::ConfigTree const& config)
+{
+    //! \ogs_file_param{material__fluid__density__type}
+    config.checkConfigParameter("type", "ConcentrationAndPressureDependent");
+
+    const double reference_density =
+        //! \ogs_file_param{material__fluid__density__ConcentrationAndPressureDependent__reference_density}
+        config.getConfigParameter<double>("reference_density");
+    const double reference_concentration =
+        //! \ogs_file_param{material__fluid__density__ConcentrationAndPressureDependent__reference_concentration}
+        config.getConfigParameter<double>("reference_concentration");
+    const double fluid_density_concentration_difference_ratio =
+        //! \ogs_file_param{material__fluid__density__ConcentrationAndPressureDependent__fluid_density_concentration_difference_ratio}
+        config.getConfigParameter<double>(
+            "fluid_density_concentration_difference_ratio");
+    const double reference_pressure =
+        //! \ogs_file_param{material__fluid__density__ConcentrationAndPressureDependent__reference_pressure}
+        config.getConfigParameter<double>("reference_pressure");
+    const double fluid_density_pressure_difference_ratio =
+        //! \ogs_file_param{material__fluid__density__ConcentrationAndPressureDependent__fluid_density_pressure_difference_ratio}
+        config.getConfigParameter<double>(
+            "fluid_density_pressure_difference_ratio");
+    return std::make_unique<LinearConcentrationAndPressureDependentDensity>(
+        reference_density,
+        reference_concentration,
+        fluid_density_concentration_difference_ratio,
+        reference_pressure,
+        fluid_density_pressure_difference_ratio);
+}
 
 std::unique_ptr<FluidProperty> createFluidDensityModel(
     BaseLib::ConfigTree const& config)
@@ -112,6 +144,8 @@ std::unique_ptr<FluidProperty> createFluidDensityModel(
         return createLinearTemperatureDependentDensity(config);
     if (type == "ConcentrationDependent")
         return createLinearConcentrationDependentDensity(config);
+    if (type == "ConcentrationAndPressureDependent")
+        return createLinearConcentrationAndPressureDependentDensity(config);
     if (type == "IdealGasLaw")
     {
         //! \ogs_file_param{material__fluid__density__type}

MaterialLib/Fluid/Density/LinearConcentrationAndPressureDependentDensity.h

@@ -0,0 +1,111 @@
+/**
+ * @copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <string>
+
+#include "BaseLib/ConfigTree.h"
+
+#include "MaterialLib/Fluid/FluidProperty.h"
+
+namespace MaterialLib
+{
+namespace Fluid
+{
+/// Linear concentration dependent density model.
+/// \f[ \varrho = \varrho_{\text{ref}}
+/// (1 + \bar \alpha (C - C_{\text{ref}}) + \bar \beta (p - p_{\text{ref}}) )
+/// \f] where
+/// - \f$ \varrho_{\text{ref}}\f$ is the reference density
+/// - \f$ \bar \alpha\f$ is the fluid density concentration difference ratio
+/// - \f$ C_{\text{ref}}\f$ is the reference concentration
+/// - \f$ \bar \beta\f$ is the fluid density pressure difference ratio
+/// - \f$ p_{\text{ref}}\f$ is the reference pressure
+class LinearConcentrationAndPressureDependentDensity final
+    : public FluidProperty
+{
+public:
+    /**
+     * @param reference_density  \f$\rho_0\f$
+     * @param reference_concentration \f$C_0\f$
+     * @param fluid_density_concentration_difference_ratio  \f$ \bar \alpha \f$
+     * in reference
+     * @param reference_pressure \f$p_0\f$
+     * @param fluid_density_pressure_difference_ratio  \f$ \bar \beta \f$ in
+     * reference Coupled groundwater flow and transport: 2. Thermohaline and 3D
+     * convection systems
+     */
+    explicit LinearConcentrationAndPressureDependentDensity(
+        const double reference_density, double reference_concentration,
+        const double fluid_density_concentration_difference_ratio,
+        double reference_pressure,
+        const double fluid_density_pressure_difference_ratio)
+        : _reference_density(reference_density),
+          _reference_concentration(reference_concentration),
+          _fluid_density_concentration_difference_ratio(
+              fluid_density_concentration_difference_ratio),
+          _reference_pressure(reference_pressure),
+          _fluid_density_pressure_difference_ratio(
+              fluid_density_pressure_difference_ratio)
+    {
+    }
+
+    /// Get model name.
+    std::string getName() const override
+    {
+        return "Linear concentration and pressure dependent density";
+    }
+
+    /// Get density value.
+    /// \param var_vals Variable values in an array. The order of its elements
+    ///                 is given in enum class PropertyVariableType.
+    double getValue(const ArrayType& var_vals) const override
+    {
+        const double C = var_vals[static_cast<int>(PropertyVariableType::C)];
+        const double p = var_vals[static_cast<int>(PropertyVariableType::p)];
+        return _reference_density *
+               (1 +
+                _fluid_density_concentration_difference_ratio *
+                    (C - _reference_concentration) +
+                _fluid_density_pressure_difference_ratio *
+                    (p - _reference_pressure));
+    }
+
+    /// Get the partial differential of the density with respect to
+    /// concentration or pressure.
+    double getdValue(const ArrayType& /*var_vals*/,
+                     const PropertyVariableType var) const override
+    {
+        if (var == PropertyVariableType::C)
+        {
+            return _reference_density *
+                   _fluid_density_concentration_difference_ratio;
+        }
+        else if (var == PropertyVariableType::p)
+        {
+            return _reference_density *
+                   _fluid_density_pressure_difference_ratio;
+        }
+        else
+        {
+            return 0;
+        }
+    }
+
+private:
+    const double _reference_density;
+    const double _reference_concentration;
+    const double _fluid_density_concentration_difference_ratio;
+    const double _reference_pressure;
+    const double _fluid_density_pressure_difference_ratio;
+};
+
+}  // namespace Fluid
+}  // namespace MaterialLib

ProcessLib/ComponentTransport/ComponentTransportFEM.h

@@ -12,7 +12,6 @@
 #include <Eigen/Dense>
 #include <vector>
 
-
 #include "ComponentTransportProcessData.h"
 #include "MaterialLib/Fluid/FluidProperties/FluidProperties.h"
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
@@ -141,7 +140,8 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface
         auto const num_nodes = ShapeFunction::NPOINTS;
         auto p_nodal_values =
             Eigen::Map<const NodalVectorType>(&local_x[num_nodes], num_nodes);
-
+        auto C_nodal_values =
+            Eigen::Map<const NodalVectorType>(&local_x[0], num_nodes);
         auto const& b = _process_data.specific_body_force;
 
         MaterialLib::Fluid::FluidProperty::ArrayType vars;
@@ -151,21 +151,18 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface
 
         auto KCC = local_K.template block<num_nodes, num_nodes>(0, 0);
         auto MCC = local_M.template block<num_nodes, num_nodes>(0, 0);
+        auto MCp = local_M.template block<num_nodes, num_nodes>(0, num_nodes);
         auto Kpp =
             local_K.template block<num_nodes, num_nodes>(num_nodes, num_nodes);
         auto Mpp =
             local_M.template block<num_nodes, num_nodes>(num_nodes, num_nodes);
+        auto MpC = local_M.template block<num_nodes, num_nodes>(num_nodes, 0);
         auto Bp = local_b.template block<num_nodes, 1>(num_nodes, 0);
 
         for (std::size_t ip(0); ip < n_integration_points; ++ip)
         {
             pos.setIntegrationPoint(ip);
 
-            // \todo the argument to getValue() has to be changed for non
-            // constant storage model
-            auto const specific_storage =
-                _process_data.porous_media_properties.getSpecificStorage(t, pos)
-                    .getValue(0.0);
 
             auto const& ip_data = _ip_data[ip];
             auto const& N = ip_data.N;
@@ -177,7 +174,6 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface
             // Order matters: First C, then p!
             NumLib::shapeFunctionInterpolate(local_x, N, C_int_pt, p_int_pt);
 
-            // \todo the first argument has to be changed for non constant
             // porosity model
             auto const porosity =
                 _process_data.porous_media_properties.getPorosity(t, pos)
@@ -210,13 +206,21 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface
                 MaterialLib::Fluid::FluidPropertyType::Viscosity, vars);
 
             GlobalDimMatrixType const K_over_mu = K / mu;
-
             GlobalDimVectorType const velocity =
                 _process_data.has_gravity
                     ? GlobalDimVectorType(-K_over_mu *
                                           (dNdx * p_nodal_values - density * b))
                     : GlobalDimVectorType(-K_over_mu * dNdx * p_nodal_values);
 
+            const double drho_dp = _process_data.fluid_properties->getdValue(
+                MaterialLib::Fluid::FluidPropertyType::Density,
+                vars,
+                MaterialLib::Fluid::PropertyVariableType::p);
+
+            const double drho_dC = _process_data.fluid_properties->getdValue(
+                MaterialLib::Fluid::FluidPropertyType::Density,
+                vars,
+                MaterialLib::Fluid::PropertyVariableType::C);
             double const velocity_magnitude = velocity.norm();
             GlobalDimMatrixType const hydrodynamic_dispersion =
                 velocity_magnitude != 0.0
@@ -236,20 +240,25 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface
                           I);
 
             // matrix assembly
+            MCp.noalias() += w * N.transpose() * N * C_nodal_values *
+                             retardation_factor * porosity * drho_dp * N;
+            MCC.noalias() += w * N.transpose() * retardation_factor * porosity *
+                                 density * N +
+                             w * N.transpose() * retardation_factor * porosity *
+                                 N * C_nodal_values * drho_dC * N;
+
             KCC.noalias() +=
-                (dNdx.transpose() * hydrodynamic_dispersion * dNdx +
-                 N.transpose() * velocity.transpose() * dNdx +
+                (-dNdx.transpose() * velocity * density * N +
+                 dNdx.transpose() * density * hydrodynamic_dispersion * dNdx +
                  N.transpose() * decay_rate * porosity * retardation_factor *
-                     N) *
+                     density * N) *
                 w;
-            MCC.noalias() +=
-                w * N.transpose() * porosity * retardation_factor * N;
-            Kpp.noalias() += w * dNdx.transpose() * K_over_mu * dNdx;
-            Mpp.noalias() += w * N.transpose() * specific_storage * N;
+
+            Mpp.noalias() += w * N.transpose() * porosity * drho_dp * N;
+            MpC.noalias() += w * N.transpose() * porosity * drho_dC * N;
+            Kpp.noalias() += w * dNdx.transpose() * density * K_over_mu * dNdx;
             if (_process_data.has_gravity)
-                Bp += w * density * dNdx.transpose() * K_over_mu * b;
-            /* with Oberbeck-Boussing assumption density difference only exists
-             * in buoyancy effects */
+                Bp += w * density * density * dNdx.transpose() * K_over_mu * b;
         }
     }
 
@@ -275,9 +284,9 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface
         pos.setElementID(_element.getID());
 
         MaterialLib::Fluid::FluidProperty::ArrayType vars;
-
-        auto const p_nodal_values = Eigen::Map<const NodalVectorType>(
-            &local_x[ShapeFunction::NPOINTS], ShapeFunction::NPOINTS);
+        auto const num_nodes = ShapeFunction::NPOINTS;
+        auto const p_nodal_values =
+            Eigen::Map<const NodalVectorType>(&local_x[num_nodes], num_nodes);
 
         for (unsigned ip = 0; ip < n_integration_points; ++ip)
         {

ProcessLib/ComponentTransport/ComponentTransportProcess.h

@@ -25,26 +25,35 @@ namespace ComponentTransport
  *
  * The flow process is described by
  * \f[
- * S \frac{\partial p}{\partial t}
- *     - \nabla \cdot \left[\frac{\kappa}{\mu(C)} \nabla \left( p + \rho g z \right)\right]
- *     - Q_p = 0,
+ * \phi \frac{\partial \rho}{\partial p} \frac{\partial p}{\partial t}
+ *     + \phi \frac{\partial \rho}{\partial C} \frac{\partial C}{\partial t}
+ *     - \nabla \cdot \left[\frac{\kappa}{\mu(C)} \rho \nabla \left( p + \rho g
+ * z \right)\right]
+ *     + Q_p = 0,
  * \f]
- * where \f$S\f$ is the storage, \f$p\f$ is the pressure,
+ * where the storage \f$S\f$ has been substituted by
+ *      \f$\phi \frac{\partial \rho}{\partial p}\f$,
+ * \f$\phi\f$ is the porosity,
+ * \f$C\f$ is the concentration,
+ * \f$p\f$ is the pressure,
  * \f$\kappa\f$ is permeability,
  * \f$\mu\f$ is viscosity of the fluid,
  * \f$\rho\f$ is the density of the fluid, and
  * \f$g\f$ is the gravitational acceleration.
  *
  * The mass transport process is described by
  * \f[
- * \phi R \frac{\partial C}{\partial t}
-    + \nabla \cdot \left(\vec{q} C - D \nabla C \right)
-        + \phi R \vartheta C - Q_C = 0
+ * \phi R C \frac{\partial \rho}{\partial p} \frac{\partial p}{\partial t}
+ *     + \phi R \left(\rho + C \frac{\partial \rho}{\partial C}\right)
+ * \frac{\partial C}{\partial t}
+ *     - \nabla \cdot \left[\frac{\kappa}{\mu(C)} \rho C \nabla \left( p + \rho
+ * g z \right)
+ *          + \rho D \nabla C\right]
+ *     + Q_C + R \vartheta \phi \rho C = 0,
+ *
  * \f]
- * where \f$\phi\f$ is the porosity,
- * \f$R\f$ is the retardation factor,
- * \f$C\f$ is the concentration,
- * \f$\vec{q} = \frac{\kappa}{\mu(C)} \nabla \left( p + \rho g z \right)\f$
+ * where \f$R\f$ is the retardation factor,
+ * \f$\vec{q} = -\frac{\kappa}{\mu(C)} \nabla \left( p + \rho g z \right)\f$
  *      is the Darcy velocity,
  * \f$D\f$ is the hydrodynamic dispersion tensor,
  * \f$\vartheta\f$ is the decay rate.
@@ -70,9 +79,9 @@ namespace ComponentTransport
  * flow couples the H process to the C process.
  *
  * \note At the moment there is not any coupling by source or sink terms, i.e.,
- * the coupling is implemented only by density changes due to concentration
- * changes in the buoyance term in the groundwater flow. This coupling schema is
- * referred to as the Boussinesq approximation.
+ * the coupling is implemented only by density and viscosity changes due to
+ * concentration changes as well as by the temporal derivatives of each
+ * variable.
  * */
 class ComponentTransportProcess final : public Process
 {