Merge branch 'iapws-thermal-conductivity' into 'master'

[MPL] Adding IAPWS model of thermal conductivity of water to MPL

See merge request !4445

Documentation/ProjectFile/properties/property/t_WaterThermalConductivityIAPWS.md

+\copydoc MaterialPropertyLib::WaterThermalConductivityIAPWS

MaterialLib/MPL/CreateProperty.cpp

@@ -300,16 +300,21 @@ std::unique_ptr<MaterialPropertyLib::Property> createProperty(
     {
         return createVolumeFractionAverage(config);
     }
+
     if (boost::iequals(property_type, "WaterViscosityIAPWS"))
     {
         return createWaterViscosityIAPWS(config);
     }
+
     if (boost::iequals(property_type, "LiquidViscosityVogels"))
     {
         return createLiquidViscosityVogels(config);
     }
 
-
+    if (boost::iequals(property_type, "WaterThermalConductivityIAPWS"))
+    {
+        return createWaterThermalConductivityIAPWS(config);
+    }
 
     // If none of the above property types are found, OGS throws an error.
     OGS_FATAL("The specified component property type '{:s}' was not recognized",

MaterialLib/MPL/Properties/CreateProperties.h

@@ -65,3 +65,4 @@
 #include "VapourDiffusion/CreateVapourDiffusionPMQ.h"
 #include "Viscosity/CreateWaterViscosityIAPWS.h"
 #include "Viscosity/CreateLiquidViscosityVogels.h"
+#include "ThermalConductivity/CreateWaterThermalConductivityIAPWS.h"

MaterialLib/MPL/Properties/ThermalConductivity/CreateWaterThermalConductivityIAPWS.cpp

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ * Created on March 4, 2021, 4:38 PM
+ */
+
+#include "CreateWaterThermalConductivityIAPWS.h"
+
+#include "BaseLib/ConfigTree.h"
+#include "MaterialLib/MPL/Property.h"
+#include "WaterThermalConductivityIAPWS.h"
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<Property> createWaterThermalConductivityIAPWS(
+    BaseLib::ConfigTree const& config)
+{
+    //! \ogs_file_param{properties__property__type}
+    config.checkConfigParameter("type", "WaterThermalConductivityIAPWS");
+    DBUG("Create WaterThermalConductivityIAPWS phase property");
+
+    // Second access for storage.
+    //! \ogs_file_param{properties__property__name}
+    auto property_name = config.peekConfigParameter<std::string>("name");
+
+    //! \ogs_file_param_special{properties__property__WaterThermalConductivityIAPWS}
+    return std::make_unique<WaterThermalConductivityIAPWS>(std::move(property_name));
+}
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/ThermalConductivity/CreateWaterThermalConductivityIAPWS.h

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ * Created on March 4, 2021, 4:38 PM
+ */
+
+#pragma once
+
+#include <memory>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace MaterialPropertyLib
+{
+class Property;
+std::unique_ptr<Property> createWaterThermalConductivityIAPWS(
+    BaseLib::ConfigTree const& config);
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/ThermalConductivity/WaterThermalConductivityIAPWS.cpp

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ * Created on March 4, 2021, 3:05 PM
+ */
+
+#include "WaterThermalConductivityIAPWS.h"
+
+#include <cmath>
+
+#include "BaseLib/Error.h"
+#include "MaterialLib/MPL/Medium.h"
+
+namespace MaterialPropertyLib
+{
+    // Li and Lij are coefficients from Tables 1 and 2 (Daucik and Dooley, 2011)
+static const double Li[5] = {2.443221e-3, 1.323095e-2, 6.770357e-3, -3.454586e-3, 4.096266e-4};
+static const double Lij[5][6] = {
+    {1.60397357, -0.646013523, 0.111443906, 0.102997357, -0.0504123634, 0.00609859258},
+    {2.33771842, -2.78843778, 1.53616167, -0.463045512, 0.0832827019, -0.00719201245},
+    {2.19650529, -4.54580785, 3.55777244, -1.40944978, 0.275418278, -0.0205938816},
+    {-1.21051378, 1.60812989, -0.621178141, 0.0716373224, 0, 0},
+    {-2.7203370, 4.57586331, -3.18369245, 1.1168348, -0.19268305, 0.012913842}};
+
+static double computeBarLambda0Factor(const double barT);
+
+static std::array<double, 5> computeSeriesFactorTForLambda1(const double barT);
+static std::array<double, 6> computeSeriesFactorRhoForLambda1(const double bar_rho);
+static double computeBarLambda1Factor(
+    const std::array<double, 5>& series_factorT,
+    const std::array<double, 6>& series_factorRho);
+
+static double computedBarLambda_dbarT(const double barT, double bar_rho);
+static double computedBarLambda_dbarRho(const double barT, double bar_rho);
+
+
+PropertyDataType WaterThermalConductivityIAPWS::value(
+    VariableArray const& variable_array,
+    ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
+    double const /*dt*/) const
+{
+    const double bar_T = variable_array.temperature / ref_T_;
+    const double bar_rho = variable_array.density / ref_rho_;
+
+    const double lambda0 = std::sqrt(bar_T) / computeBarLambda0Factor(bar_T);
+
+    const auto& series_factorT = computeSeriesFactorTForLambda1(bar_T);
+    const auto& series_factorRho = computeSeriesFactorRhoForLambda1(bar_rho);
+    const double lambda1 = std::exp(
+        bar_rho * computeBarLambda1Factor(series_factorT, series_factorRho));
+
+    return lambda0 * lambda1 * ref_lambda_;
+}
+
+PropertyDataType WaterThermalConductivityIAPWS::dValue(
+    VariableArray const& variable_array, Variable const variable,
+    ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
+    double const /*dt*/) const
+{
+    const double bar_T = variable_array.temperature / ref_T_;
+    const double bar_rho = variable_array.density / ref_rho_;
+
+    switch (variable)
+    {
+        case Variable::temperature:
+            return ref_lambda_ * computedBarLambda_dbarT(bar_T, bar_rho) / ref_T_;
+        case Variable::density:
+            return ref_lambda_ * computedBarLambda_dbarRho(bar_T, bar_rho) / ref_rho_;
+        default:
+            OGS_FATAL(
+            "WaterThermalConductivityIAPWS::dValue is implemented for "
+            "derivatives with respect to temperature and liquid density only.");
+    }
+}
+
+double computeBarLambda0Factor(const double barT)
+{
+    double sum_val = 0.;
+    double barT_i = 1.;
+    for (double value : Li)
+    {
+        sum_val += (value / barT_i);
+        barT_i *= barT;
+    }
+    return sum_val;
+}
+
+std::array<double, 5> computeSeriesFactorTForLambda1(const double barT)
+{
+    std::array<double, 5> series_factorT;
+    series_factorT[0] = 1.;
+    const double barT_fac = 1 / barT - 1.0;
+    for (int i = 1; i < 5; i++)
+    {
+        series_factorT[i] = series_factorT[i - 1] * barT_fac;
+    }
+
+    return series_factorT;
+}
+
+std::array<double, 6> computeSeriesFactorRhoForLambda1(const double bar_rho)
+{
+    std::array<double, 6> series_factorRho;
+    series_factorRho[0] = 1.;
+    for (int i = 1; i < 6; i++)
+    {
+        series_factorRho[i] = series_factorRho[i - 1] * (bar_rho - 1.0);
+    }
+
+    return series_factorRho;
+}
+
+double computeBarLambda1Factor(const std::array<double, 5>& series_factorT,
+                           const std::array<double, 6>& series_factorRho)
+{
+    double sum_val = 0.;
+    for (int i = 0; i < 5; i++)
+    {
+        double sum_val_j = 0;
+        for (int j = 0; j < 6; j++)
+        {
+            sum_val_j += Lij[i][j] * series_factorRho[j];
+        }
+        sum_val += series_factorT[i] * sum_val_j;
+    }
+
+    return sum_val;
+}
+
+double computedBarLambda_dbarT(const double barT, double bar_rho)
+{
+    const double lambda0_factor = computeBarLambda0Factor(barT);
+    const double sqrt_barT = std::sqrt(barT);
+
+    double dlambda0_factor_dbarT = 0.0;
+    double barT_i = barT * barT;
+    for (int i = 1; i < 5; i++)
+    {
+        dlambda0_factor_dbarT -= static_cast<double>(i) * (Li[i] / barT_i);
+        barT_i *= barT;
+    }
+
+    const double dbar_lambda0_dbarT =
+        0.5 / (lambda0_factor * sqrt_barT) -
+        sqrt_barT * dlambda0_factor_dbarT / (lambda0_factor * lambda0_factor);
+
+    const auto& series_factorT = computeSeriesFactorTForLambda1(barT);
+    const auto& series_factorRho = computeSeriesFactorRhoForLambda1(bar_rho);
+
+    double dlambda1_factor_dbarT = 0.0;
+    for (int i = 1; i < 5; i++)
+    {
+        double sum_val_j = 0;
+        for (int j = 0; j < 6; j++)
+        {
+            sum_val_j += Lij[i][j] * series_factorRho[j];
+        }
+        dlambda1_factor_dbarT -= static_cast<double>(i) * series_factorT[i - 1] *
+                             sum_val_j / (barT * barT);
+    }
+
+    const double lambda1_factor =
+        computeBarLambda1Factor(series_factorT, series_factorRho);
+    const double dbar_lambda1_dbarT =
+        bar_rho * std::exp(bar_rho * lambda1_factor) * dlambda1_factor_dbarT;
+
+    return dbar_lambda0_dbarT * std::exp(bar_rho * lambda1_factor) +
+           dbar_lambda1_dbarT * sqrt_barT / lambda0_factor;
+}
+
+double computedBarLambda_dbarRho(const double barT, double bar_rho)
+{
+    const auto& series_factorT = computeSeriesFactorTForLambda1(barT);
+    const auto& series_factorRho = computeSeriesFactorRhoForLambda1(bar_rho);
+
+    double dlambda1_factor_dbar_rho = 0.0;
+    for (int i = 0; i < 5; i++)
+    {
+        double sum_val_j = 0;
+        for (int j = 1; j < 6; j++)
+        {
+            sum_val_j +=
+                static_cast<double>(j) * Lij[i][j] * series_factorRho[j - 1];
+        }
+        dlambda1_factor_dbar_rho += series_factorT[i] * sum_val_j;
+    }
+
+    const double lambda0 = std::sqrt(barT) / computeBarLambda0Factor(barT);
+
+    const double lambda1_factor =
+        computeBarLambda1Factor(series_factorT, series_factorRho);
+    return lambda0 * std::exp(bar_rho * lambda1_factor) *
+           (lambda1_factor + bar_rho * dlambda1_factor_dbar_rho);
+}
+
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/ThermalConductivity/WaterThermalConductivityIAPWS.h

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ * Created on March 4, 2021, 3:05 PM
+ */
+
+#pragma once
+
+#include "MaterialLib/MPL/Property.h"
+
+namespace MaterialPropertyLib
+{
+class Phase;
+
+/**
+ * \brief A class for thermal conductivity model that is defined by
+ *        The International Association for the Properties of Water and Steam
+ *        <a href="http://www.iapws.org/relguide/ThCond.pdf">IAPWS</a>
+ *        (File accessed at 13.01.2023) - (Daucik and Dooley, 2011)
+ *
+ *        With the definition, the thermal conductivity is a function of temperature and
+ *        water density
+ *
+ *  \attention The critical enhancement, \f$\bar{\lambda}_2\f$, is not considered.
+ *              For information on region of significance and the significance,
+ *              please see Figure 2 from the document linked in the upper paragraph.
+ */
+class WaterThermalConductivityIAPWS final : public Property
+{
+public:
+    explicit WaterThermalConductivityIAPWS(std::string name) { name_ = std::move(name); }
+    void checkScale() const override
+    {
+        if (!std::holds_alternative<Phase*>(scale_))
+        {
+            OGS_FATAL(
+                "The property 'WaterThermalConductivityIAPWS' is "
+                "implemented on the 'Phase' scale only.");
+        }
+    }
+
+    /// \return The water thermal conductivity.
+    PropertyDataType value(VariableArray const& variable_array,
+                           ParameterLib::SpatialPosition const& pos,
+                           double const t, double const dt) const override;
+    /// \return The derivative of water thermal conductivity with respect to
+    /// temperature or phase (water) pressure.
+    PropertyDataType dValue(VariableArray const& variable_array,
+                            Variable const variable,
+                            ParameterLib::SpatialPosition const& pos,
+                            double const t, double const dt) const override;
+private:
+    static constexpr double ref_T_ = 647.096;  ///< reference temperature in K
+    static constexpr double ref_rho_ = 322.0;  ///< reference density in `kg/m^3`
+    static constexpr double ref_lambda_ = 1.0e-3;  ///< reference thermal conductivity in `W.K^-1.m^-1`
+
+};
+}  // namespace MaterialPropertyLib

Tests/MaterialLib/TestMPLWaterThermalConductivity.cpp

+/**
+ *  \file
+ *  \brief Test thermal conductivity of water IAPWS97 model
+ *
+ *  \copyright
+ *   Copyright (c) 2012-2023, 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 <memory>
+
+#include "MaterialLib/MPL/Properties/ThermalConductivity/CreateWaterThermalConductivityIAPWS.h"
+#include "TestMPL.h"
+
+TEST(Material, checkWaterThermalConductivityIAPWS_)
+{
+    const char xml[] =
+        "<property>"
+        "  <name>thermal_conductivity</name>"
+        "  <type>WaterThermalConductivityIAPWS</type>"
+        "</property>";
+
+    std::unique_ptr<MaterialPropertyLib::Property> const property_ptr =
+        Tests::createTestProperty(
+            xml, MaterialPropertyLib::createWaterThermalConductivityIAPWS);
+    MaterialPropertyLib::Property const& property = *property_ptr;
+
+    MaterialPropertyLib::VariableArray variable_array;
+    ParameterLib::SpatialPosition const pos;
+    double const t = std::numeric_limits<double>::quiet_NaN();
+    double const dt = std::numeric_limits<double>::quiet_NaN();
+
+    // Test data provided on http://www.iapws.org/relguide/ThCond.pdf Table 4
+    const double T[] = {298.15, 298.15, 298.15, 873.15};
+    const double rho[] = {0, 998, 1200, 0};
+
+    const double expected_lambda[] = {18.4341883e-3, 607.712868e-3, 799.038144e-3, 79.1034659e-3};
+
+    const double perturbation = 1.e-8;
+    for (int i = 0; i < 4; i++)
+    {
+        // Test lambda
+        variable_array.temperature = T[i];
+        variable_array.density = rho[i];
+        const double lambda = property.template value<double>(
+                variable_array, pos, t, dt);
+        ASSERT_NEAR(expected_lambda[i], lambda, 1.e-9);
+
+        const double dlambda_dT = property.template dValue<double>(
+                variable_array, MaterialPropertyLib::Variable::temperature, pos,
+                t, dt);
+        const double dlambda_drho = property.template dValue<double>(
+                variable_array, MaterialPropertyLib::Variable::density, pos,
+                t, dt);
+
+        // Test dlambda/dT
+        variable_array.temperature =
+            T[i] + perturbation;
+        double lambda1 = property.template value<double>(
+                variable_array, pos, t, dt);
+        ASSERT_NEAR((lambda1 - lambda) / perturbation, dlambda_dT, 8.e-6);
+
+        // Test dlambda/drho
+        variable_array.temperature = T[i];
+        variable_array.density =
+            rho[i] + perturbation;
+        lambda1 = property.template value<double>(
+                variable_array, pos, t, dt);
+
+        ASSERT_NEAR((lambda1 - lambda) / perturbation, dlambda_drho, 8.e-6);
+    }
+}