[MPL] medium effective thermal conductivity.

Documentation/ProjectFile/properties/property/EffectiveThermalConductivityPorosityMixing/c_EffectiveThermalConductivityPorosityMixing.md

+\copydoc MaterialPropertyLib::EffectiveThermalConductivityPorosityMixing

MaterialLib/MPL/CreatePhase.h

@@ -15,6 +15,7 @@
 #include <map>
 #include <memory>
 #include <optional>
+#include <string>
 #include <vector>
 
 namespace BaseLib

MaterialLib/MPL/CreateProperty.cpp

@@ -76,6 +76,12 @@ std::unique_ptr<MaterialPropertyLib::Property> createProperty(
         return createDupuitPermeability(config, parameters);
     }
 
+    if (boost::iequals(property_type, "EffectiveThermalConductivityPorosityMixing"))
+    {
+        return createEffectiveThermalConductivityPorosityMixing(
+                geometry_dimension, config, local_coordinate_system);
+    }
+
     if (boost::iequals(property_type, "IdealGasLaw"))
     {
         return createIdealGasLaw(config);

MaterialLib/MPL/Properties/CreateEffectiveThermalConductivityPorosityMixing.cpp

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ */
+
+#include "BaseLib/ConfigTree.h"
+#include "MaterialLib/MPL/Property.h"
+#include "ParameterLib/CoordinateSystem.h"
+#include "EffectiveThermalConductivityPorosityMixing.h"
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<Property> createEffectiveThermalConductivityPorosityMixing(
+        int const geometry_dimension,
+        BaseLib::ConfigTree const& config,
+        ParameterLib::CoordinateSystem const* const local_coordinate_system)
+{
+    //! \ogs_file_param{properties__property__type}
+    config.checkConfigParameter("type", "EffectiveThermalConductivityPorosityMixing");
+
+    // Second access for storage.
+    //! \ogs_file_param{properties__property__name}
+    auto property_name = config.peekConfigParameter<std::string>("name");
+
+    DBUG("Create effective thermal_conductivity property from porosity mixing {:s}.",
+         property_name);
+
+
+    if (geometry_dimension == 1)
+    {
+        return std::make_unique<
+            MaterialPropertyLib::EffectiveThermalConductivityPorosityMixing<1>>(
+            std::move(property_name),
+            local_coordinate_system);
+    }
+
+    if (geometry_dimension == 2)
+    {
+        return std::make_unique<
+            MaterialPropertyLib::EffectiveThermalConductivityPorosityMixing<2>>(
+            std::move(property_name),
+            local_coordinate_system);
+    }
+    //! \ogs_file_param_special{properties__property__EffectiveThermalConductivityPorosityMixing}
+    return std::make_unique<
+            MaterialPropertyLib::EffectiveThermalConductivityPorosityMixing<3>>(
+            std::move(property_name),
+            local_coordinate_system);
+}
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/CreateEffectiveThermalConductivityPorosityMixing.h

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, 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 <memory>
+
+namespace MaterialPropertyLib
+{
+class Property;
+
+std::unique_ptr<Property> createEffectiveThermalConductivityPorosityMixing(
+        int const geometry_dimension,
+        BaseLib::ConfigTree const& config,
+        ParameterLib::CoordinateSystem const* const local_coordinate_system);
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/CreateProperties.h

@@ -23,6 +23,7 @@
 #include "CreateConstant.h"
 #include "CreateCurve.h"
 #include "CreateDupuitPermeability.h"
+#include "CreateEffectiveThermalConductivityPorosityMixing.h"
 #include "CreateEmbeddedFracturePermeability.h"
 #include "CreateExponential.h"
 #include "CreateGasPressureDependentPermeability.h"

MaterialLib/MPL/Properties/EffectiveThermalConductivityPorosityMixing.cpp

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ */
+
+#include "EffectiveThermalConductivityPorosityMixing.h"
+
+#include "MaterialLib/MPL/Medium.h"
+#include "MaterialLib/MPL/Utils/FormEigenTensor.h"
+#include "ParameterLib/CoordinateSystem.h"
+#include "ParameterLib/Parameter.h"
+
+namespace MaterialPropertyLib
+{
+//
+// For 1D problems
+//
+template <>
+EffectiveThermalConductivityPorosityMixing<1>::EffectiveThermalConductivityPorosityMixing(
+    std::string name,
+    ParameterLib::CoordinateSystem const* const local_coordinate_system)
+    : local_coordinate_system_(local_coordinate_system)
+{
+    name_ = std::move(name);
+}
+
+template <>
+void EffectiveThermalConductivityPorosityMixing<1>::checkScale() const
+{
+    if (!std::holds_alternative<Medium*>(scale_))
+    {
+        OGS_FATAL(
+            "The property 'EffectiveThermalConductivityPorosityMixing' is "
+            "implemented on the 'medium' scale only.");
+    }
+}
+
+template <>
+PropertyDataType EffectiveThermalConductivityPorosityMixing<1>::value(
+    VariableArray const& variable_array,
+    ParameterLib::SpatialPosition const& pos, double const t,
+    double const dt) const
+{
+    auto const& medium = std::get<Medium*>(scale_);
+    auto const& liquid_phase = medium->phase("AqueousLiquid");
+    auto const& solid_phase = medium->phase("Solid");
+    auto const porosity = medium->property(
+                        MaterialPropertyLib::PropertyType::porosity)
+                    .template value<double>(variable_array, pos, t, dt);
+    auto const liquid_thermal_conductivity =
+        liquid_phase
+            .property(MaterialPropertyLib::PropertyType::thermal_conductivity)
+            .template value<double>(variable_array, pos, t, dt);
+    auto const solid_thermal_conductivity = solid_phase.property(
+                        MaterialPropertyLib::PropertyType::thermal_conductivity)
+                    .template value<double>(variable_array, pos, t, dt);
+
+    double const effective_thermal_conductivity =
+        (1.0 - porosity) * solid_thermal_conductivity +
+        porosity * liquid_thermal_conductivity;
+    return effective_thermal_conductivity;
+}
+template <>
+PropertyDataType EffectiveThermalConductivityPorosityMixing<1>::dValue(
+    VariableArray const& variable_array, Variable const variable,
+    ParameterLib::SpatialPosition const& pos, double const t,
+    double const dt) const
+{
+    auto const& medium = std::get<Medium*>(scale_);
+    auto const& liquid_phase = medium->phase("AqueousLiquid");
+    auto const& solid_phase = medium->phase("Solid");
+    auto const porosity = medium->property(
+                        MaterialPropertyLib::PropertyType::porosity)
+                    .template value<double>(variable_array, pos, t, dt);
+    auto const liquid_thermal_conductivity =
+        liquid_phase
+            .property(MaterialPropertyLib::PropertyType::thermal_conductivity)
+            .template value<double>(variable_array, pos, t, dt);
+    auto const solid_thermal_conductivity = solid_phase.property(
+                        MaterialPropertyLib::PropertyType::thermal_conductivity)
+                    .template value<double>(variable_array, pos, t, dt);
+    auto const dporosity = medium->property(
+                        MaterialPropertyLib::PropertyType::porosity)
+                    .template dValue<double>(variable_array, variable, pos, t, dt);
+    auto const dliquid_thermal_conductivity =
+        liquid_phase
+            .property(MaterialPropertyLib::PropertyType::thermal_conductivity)
+            .template dValue<double>(variable_array, variable, pos, t, dt);
+    auto const dsolid_thermal_conductivity = solid_phase.property(
+                        MaterialPropertyLib::PropertyType::thermal_conductivity)
+                    .template dValue<double>(variable_array, variable, pos, t, dt);
+
+    double const deffective_thermal_conductivity =
+        -1.0 * dporosity * solid_thermal_conductivity +
+        (1.0 - porosity) * dsolid_thermal_conductivity +
+        dporosity * liquid_thermal_conductivity +
+        porosity * dliquid_thermal_conductivity;
+    return deffective_thermal_conductivity;
+}
+//
+// For 2D and 3D problems
+//
+template <int GlobalDim>
+EffectiveThermalConductivityPorosityMixing<GlobalDim>::EffectiveThermalConductivityPorosityMixing(
+    std::string name,
+    ParameterLib::CoordinateSystem const* const local_coordinate_system)
+    : local_coordinate_system_(local_coordinate_system)
+{
+    name_ = std::move(name);
+}
+
+template <int GlobalDim>
+void EffectiveThermalConductivityPorosityMixing<GlobalDim>::checkScale() const
+{
+    if (!std::holds_alternative<Medium*>(scale_))
+    {
+        OGS_FATAL(
+            "The property 'EffectiveThermalConductivityPorosityMixing' is "
+            "implemented on the 'medium' scale only.");
+    }
+}
+
+template <int GlobalDim>
+PropertyDataType EffectiveThermalConductivityPorosityMixing<GlobalDim>::value(
+    VariableArray const& variable_array,
+    ParameterLib::SpatialPosition const& pos, double const t,
+    double const dt) const
+{
+    auto const& medium = std::get<Medium*>(scale_);
+    auto const& liquid_phase = medium->phase("AqueousLiquid");
+    auto const& solid_phase = medium->phase("Solid");
+    auto const porosity = medium->property(
+                        MaterialPropertyLib::PropertyType::porosity)
+                    .template value<double>(variable_array, pos, t, dt);
+    auto const liquid_thermal_conductivity =
+        liquid_phase
+            .property(MaterialPropertyLib::PropertyType::thermal_conductivity)
+            .template value<double>(variable_array, pos, t, dt);
+    auto solid_thermal_conductivity = formEigenTensor<GlobalDim>(
+        solid_phase
+            .property(MaterialPropertyLib::PropertyType::thermal_conductivity)
+            .value(variable_array, pos, t, dt));
+
+    // Local coordinate transformation is only applied for the case that the
+    // initial solid thermal conductivity is given with orthotropic assumption.
+    if (local_coordinate_system_ && (solid_thermal_conductivity.cols() == GlobalDim))
+    {
+        Eigen::Matrix<double, GlobalDim, GlobalDim> const e =
+            local_coordinate_system_->transformation<GlobalDim>(pos);
+
+        solid_thermal_conductivity = e.transpose() * solid_thermal_conductivity * e;
+    }
+    Eigen::Matrix<double, GlobalDim, GlobalDim> const
+        effective_thermal_conductivity =
+            (1.0 - porosity) * solid_thermal_conductivity +
+            porosity * liquid_thermal_conductivity *
+                Eigen::Matrix<double, GlobalDim, GlobalDim>::Identity();
+    return effective_thermal_conductivity;
+}
+
+template <int GlobalDim>
+PropertyDataType EffectiveThermalConductivityPorosityMixing<GlobalDim>::dValue(
+    VariableArray const& variable_array, Variable const variable,
+    ParameterLib::SpatialPosition const& pos, double const t,
+    double const dt) const
+{
+    auto const& medium = std::get<Medium*>(scale_);
+    auto const& liquid_phase = medium->phase("AqueousLiquid");
+    auto const& solid_phase = medium->phase("Solid");
+    auto const porosity = medium->property(
+                        MaterialPropertyLib::PropertyType::porosity)
+                    .template value<double>(variable_array, pos, t, dt);
+    auto const liquid_thermal_conductivity =
+            liquid_phase
+                .property(
+                    MaterialPropertyLib::PropertyType::thermal_conductivity)
+                .template value<double>(variable_array, pos, t, dt);
+    auto solid_thermal_conductivity =
+        formEigenTensor<GlobalDim>(solid_phase
+                    .property(
+                        MaterialPropertyLib::PropertyType::thermal_conductivity)
+                    .value(variable_array, pos, t, dt));
+    auto const dporosity = medium->property(
+                        MaterialPropertyLib::PropertyType::porosity)
+                    .template dValue<double>(variable_array, variable, pos, t, dt);
+    auto const dliquid_thermal_conductivity =
+            liquid_phase
+                .property(
+                    MaterialPropertyLib::PropertyType::thermal_conductivity)
+                .template dValue<double>(variable_array, variable, pos, t, dt);
+    auto dsolid_thermal_conductivity =
+        formEigenTensor<GlobalDim>(solid_phase
+                    .property(
+                        MaterialPropertyLib::PropertyType::thermal_conductivity)
+                    .dValue(variable_array, variable, pos, t, dt));
+
+    // Local coordinate transformation is only applied for the case that the
+    // initial solid thermal conductivity is given with orthotropic assumption.
+    if (local_coordinate_system_ && (solid_thermal_conductivity.cols() == GlobalDim))
+    {
+        Eigen::Matrix<double, GlobalDim, GlobalDim> const e =
+            local_coordinate_system_->transformation<GlobalDim>(pos);
+
+        solid_thermal_conductivity = e.transpose() * solid_thermal_conductivity * e;
+        dsolid_thermal_conductivity =
+            e.transpose() * dsolid_thermal_conductivity * e;
+    }
+    Eigen::Matrix<double, GlobalDim, GlobalDim> const
+        deffective_thermal_conductivity =
+            -1.0 * dporosity * solid_thermal_conductivity +
+            (1.0 - porosity) * dsolid_thermal_conductivity +
+            dporosity * liquid_thermal_conductivity *
+                Eigen::Matrix<double, GlobalDim, GlobalDim>::Identity() +
+            porosity * dliquid_thermal_conductivity *
+                Eigen::Matrix<double, GlobalDim, GlobalDim>::Identity();
+    return deffective_thermal_conductivity;
+}
+template class EffectiveThermalConductivityPorosityMixing<2>;
+template class EffectiveThermalConductivityPorosityMixing<3>;
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/EffectiveThermalConductivityPorosityMixing.h

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, 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 "MaterialLib/MPL/Property.h"
+
+namespace ParameterLib
+{
+struct CoordinateSystem;
+template <typename T>
+struct Parameter;
+}  // namespace ParameterLib
+
+namespace MaterialPropertyLib
+{
+class Medium;
+/// Porosity mixing based model for effective heat conduction
+/// \details This property is a medium property.
+/// The corresponding values are taken from the liquid/solid phase.
+template <int GlobalDim>
+class EffectiveThermalConductivityPorosityMixing final : public Property
+{
+public:
+    EffectiveThermalConductivityPorosityMixing(std::string name,
+            ParameterLib::CoordinateSystem const* const local_coordinate_system);
+
+    void checkScale() const override;
+
+    PropertyDataType value(VariableArray const& variable_array,
+                           ParameterLib::SpatialPosition const& pos,
+                           double const t,
+                           double const dt) const override;
+    PropertyDataType dValue(VariableArray const& variable_array,
+                            Variable const variable,
+                            ParameterLib::SpatialPosition const& /*pos*/,
+                            double const /*t*/,
+                            double const /*dt*/) const override;
+
+private:
+    ParameterLib::CoordinateSystem const* const local_coordinate_system_;
+};
+
+extern template class EffectiveThermalConductivityPorosityMixing<2>;
+extern template class EffectiveThermalConductivityPorosityMixing<3>;
+
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/Properties.h

@@ -22,6 +22,7 @@
 #include "Curve.h"
 #include "Density/WaterVapourDensity.h"
 #include "DupuitPermeability.h"
+#include "EffectiveThermalConductivityPorosityMixing.h"
 #include "EmbeddedFracturePermeability.h"
 #include "Enthalpy/LinearWaterVapourLatentHeat.h"
 #include "Exponential.h"

Tests/MaterialLib/TestEffectiveThermalConductivityPorosityMixing.cpp

+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+#include <gtest/gtest.h>
+
+#include "MaterialLib/MPL/Medium.h"
+#include "MaterialLib/MPL/Properties/EffectiveThermalConductivityPorosityMixing.h"
+#include "MaterialLib/MPL/Utils/FormEigenTensor.h"
+#include "ParameterLib/ConstantParameter.h"
+#include "ParameterLib/CoordinateSystem.h"
+#include "TestMPL.h"
+
+TEST(MaterialPropertyLib, EffectiveThermalConductivityPorosityMixing)
+{
+    std::string m =
+        "<medium>"
+        "<phases><phase><type>AqueousLiquid</type>"
+        "<properties>"
+        "  <property>"
+        "    <name>thermal_conductivity</name>"
+        "    <type>Constant</type>"
+        "    <value>0.123</value>"
+        "  </property> "
+        "</properties>"
+        "</phase>"
+        "<phase><type>Solid</type>"
+        "<properties>"
+        "  <property>"
+        "    <name>thermal_conductivity</name>"
+        "    <type>Constant</type>"
+        "    <value>0.923</value>"
+        "  </property> "
+        "</properties>"
+        "</phase></phases>"
+        "<properties>"
+        "  <property>"
+        "    <name>porosity</name>"
+        "    <type>Constant</type>"
+        "    <value>0.12</value>"
+        "  </property> "
+        "  <property>"
+        "    <name>thermal_conductivity</name>"
+        "    <type>EffectiveThermalConductivityPorosityMixing</type>"
+        "  </property> "
+        "</properties>"
+        "</medium>";
+
+    auto const& medium = Tests::createTestMaterial(m, 3);
+
+    MaterialPropertyLib::VariableArray variable_array;
+    ParameterLib::SpatialPosition const pos;
+    double const time = std::numeric_limits<double>::quiet_NaN();
+    double const dt = std::numeric_limits<double>::quiet_NaN();
+
+    auto const eff_th_cond = MaterialPropertyLib::formEigenTensor<3>(
+            medium->property(MaterialPropertyLib::PropertyType::thermal_conductivity)
+                .value(variable_array, pos, time, dt));
+    auto const deff_th_cond = MaterialPropertyLib::formEigenTensor<3>(
+        medium->property(MaterialPropertyLib::PropertyType::thermal_conductivity)
+                .dValue(variable_array, MaterialPropertyLib::Variable::phase_pressure, pos, time, dt));
+    ASSERT_NEAR(eff_th_cond.trace(), 2.481, 1.e-10);
+    ASSERT_NEAR(deff_th_cond.trace(), 0.0, 1.e-10);
+}
+
+TEST(MaterialPropertyLib, EffectiveThermalConductivityPorosityMixingRot90deg)
+{
+    ParameterLib::ConstantParameter<double> const a{"a", {1., 0., 0.}};
+    ParameterLib::ConstantParameter<double> const b{"a", {0., 1., 0.}};
+    ParameterLib::ConstantParameter<double> const c{"a", {0., 0., 1.}};
+    ParameterLib::CoordinateSystem const coordinate_system{b, c, a};
+
+    std::string m =
+        "<medium>"
+        "<phases><phase><type>AqueousLiquid</type>"
+        "<properties>"
+        "  <property>"
+        "    <name>thermal_conductivity</name>"
+        "    <type>Constant</type>"
+        "    <value>0.0</value>"
+        "  </property> "
+        "</properties>"
+        "</phase>"
+        "<phase><type>Solid</type>"
+        "<properties>"
+        "  <property>"
+        "    <name>thermal_conductivity</name>"
+        "    <type>Constant</type>"
+        "    <value>0.923 0.531 0.89</value>"
+        "  </property> "
+        "</properties>"
+        "</phase></phases>"
+        "<properties>"
+        "  <property>"
+        "    <name>porosity</name>"
+        "    <type>Constant</type>"
+        "    <value>0.12</value>"
+        "  </property> "
+        "  <property>"
+        "    <name>thermal_conductivity</name>"
+        "    <type>EffectiveThermalConductivityPorosityMixing</type>"
+        "  </property> "
+        "</properties>"
+        "</medium>";
+
+    auto const& medium = Tests::createTestMaterial(m, 3, &coordinate_system);
+
+    MaterialPropertyLib::VariableArray variable_array;
+    ParameterLib::SpatialPosition const pos;
+    double const time = std::numeric_limits<double>::quiet_NaN();
+    double const dt = std::numeric_limits<double>::quiet_NaN();
+
+    auto const eff_th_cond = MaterialPropertyLib::formEigenTensor<3>(
+        medium
+            ->property(MaterialPropertyLib::PropertyType::thermal_conductivity)
+            .value(variable_array, pos, time, dt));
+    auto const deff_th_cond = MaterialPropertyLib::formEigenTensor<3>(
+        medium
+            ->property(MaterialPropertyLib::PropertyType::thermal_conductivity)
+            .dValue(variable_array,
+                    MaterialPropertyLib::Variable::phase_pressure,
+                    pos,
+                    time,
+                    dt));
+    ASSERT_NEAR(eff_th_cond(0, 0), 0.467280, 1.e-10);
+    ASSERT_NEAR(eff_th_cond(1, 1), 0.7832, 1.e-10);
+    ASSERT_NEAR(eff_th_cond(2, 2), 0.81224, 1.e-10);
+    ASSERT_NEAR(deff_th_cond.trace(), 0.0, 1.e-10);
+}

Tests/MaterialLib/TestMPL.cpp

@@ -17,13 +17,15 @@
 #include "BaseLib/ConfigTree.h"
 #include "MaterialLib/MPL/CreateMedium.h"
 #include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h"
+#include "ParameterLib/CoordinateSystem.h"
 #include "ParameterLib/Parameter.h"
 #include "Tests/TestTools.h"
 
 namespace Tests
 {
-std::unique_ptr<MPL::Medium> createTestMaterial(std::string const& xml,
-                                                int const geometry_dimension)
+std::unique_ptr<MPL::Medium> createTestMaterial(
+    std::string const& xml, int const geometry_dimension,
+    ParameterLib::CoordinateSystem const* const local_coordinate_system)
 {
     auto const ptree = Tests::readXml(xml.c_str());
     BaseLib::ConfigTree conf(ptree, "", BaseLib::ConfigTree::onerror,
@@ -34,8 +36,8 @@ std::unique_ptr<MPL::Medium> createTestMaterial(std::string const& xml,
              std::unique_ptr<MathLib::PiecewiseLinearInterpolation>>
         curves;
 
-    return MPL::createMedium(geometry_dimension, config, parameters, nullptr,
-                             curves);
+    return MPL::createMedium(geometry_dimension, config, parameters,
+                             local_coordinate_system, curves);
 }
 
 std::unique_ptr<MaterialPropertyLib::Property> createTestProperty(

Tests/MaterialLib/TestMPL.h

@@ -20,6 +20,11 @@
 
 namespace MPL = MaterialPropertyLib;
 
+namespace ParameterLib
+{
+struct CoordinateSystem;
+}
+
 namespace BaseLib
 {
 class ConfigTree;
@@ -33,7 +38,9 @@ class Property;
 namespace Tests
 {
 std::unique_ptr<MPL::Medium> createTestMaterial(
-    std::string const& xml, int const geometry_dimension = 1);
+    std::string const& xml,
+    int const geometry_dimension = 1,
+    ParameterLib::CoordinateSystem const* const local_coordinate_system = nullptr);
 
 std::unique_ptr<MaterialPropertyLib::Property> createTestProperty(
     const char xml[],