Merge branch 'MPL_vapour_diffusion' into 'master'

Added vapour diffusion type and the FEBEX type vapour_diffusion property

See merge request !3491

Documentation/ProjectFile/properties/property/VapourDiffusionFEBEX/c_VapourDiffusionFEBEX.md

+\copydoc MaterialPropertyLib::VapourDiffusionFEBEX

Documentation/ProjectFile/properties/property/VapourDiffusionFEBEX/t_tortuosity.md

+The tortuosity.

Documentation/bibliography/other.bib

@@ -228,3 +228,47 @@ URL = {https://doi.org/10.1680/geot.2008.58.3.157}
   year={1974},
   organization={Society of Petroleum Engineers}
 }
+
+@article{moldrup1997modeling,
+	title={Modeling diffusion and reaction in soils: VII. Predicting gas and ion diffusivity in undisturbed and sieved soils},
+	author={Moldrup, P and Olesen, Torben and Rolston, DE and Yamaguchi, T},
+	journal={Soil Science},
+	volume={162},
+	number={9},
+	pages={632--640},
+	year={1997},
+	publisher={LWW}
+}
+
+@article{moldrup2000predicting,
+	title={Predicting the gas diffusion coefficient in undisturbed soil from soil water characteristics},
+	author={Moldrup, P and Olesen, Torben and Schj{\o}nning, Per and Yamaguchi, T and Rolston, DE},
+	journal={Soil Science Society of America Journal},
+	volume={64},
+	number={1},
+	pages={94--100},
+	year={2000},
+	publisher={Wiley Online Library}
+}
+@article{chau2005simulation,
+	title={Simulation of gaseous diffusion in partially saturated porous media under variable gravity with lattice Boltzmann methods},
+	author={Chau, Jessica Furrer and Or, Dani and Sukop, Michael C},
+	journal={Water resources research},
+	volume={41},
+	number={8},
+	year={2005},
+	publisher={Wiley Online Library}
+}
+
+
+@inproceedings{Rutquist2007TaskA1,
+	author = {Rutquist, J.},
+	editor = {Jing L. and Nguyen, S.},
+	title = {{SKI/LBNL Modelling of Task A-1 using ROCMAS}},
+	booktitle = {Influence of near field
+	coupled {THM} phenomena on the performance of a spent fuel repository},
+	date = {2007},
+	volumes = {DECOVALEX-THMC Project, ISSN 1104-1374},
+	series = {SKI Report 2007:07},
+	pages = {61-72},
+}

MaterialLib/CMakeLists.txt

@@ -20,6 +20,7 @@ append_source_files(SOURCES MPL/Properties/CapillaryPressureSaturation)
 append_source_files(SOURCES MPL/Properties/RelativePermeability)
 append_source_files(SOURCES MPL/Properties/SwellingStress)
 append_source_files(SOURCES MPL/Properties/ThermalConductivity)
+append_source_files(SOURCES MPL/Properties/VapourDiffusion)
 append_source_files(SOURCES MPL/Components)
 append_source_files(SOURCES MPL/Utils)
 

MaterialLib/MPL/CreateProperty.cpp

@@ -221,6 +221,11 @@ std::unique_ptr<MaterialPropertyLib::Property> createProperty(
             geometry_dimension, config, parameters, local_coordinate_system);
     }
 
+    if (boost::iequals(property_type, "VapourDiffusionFEBEX"))
+    {
+        return createVapourDiffusionFEBEX(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",
               property_type);

MaterialLib/MPL/Properties/CreateProperties.h

@@ -46,3 +46,4 @@
 #include "RelativePermeability/CreateRelPermVanGenuchten.h"
 #include "SwellingStress/CreateLinearSaturationSwellingStress.h"
 #include "ThermalConductivity/CreateSoilThermalConductivitySomerton.h"
+#include "VapourDiffusion/CreateVapourDiffusionFEBEX.h"

MaterialLib/MPL/Properties/Properties.h

@@ -39,3 +39,4 @@
 #include "SaturationDependentSwelling.h"
 #include "ThermalConductivity/SoilThermalConductivitySomerton.h"
 #include "TransportPorosityFromMassBalance.h"
+#include "VapourDiffusion/VapourDiffusionFEBEX.h"

MaterialLib/MPL/Properties/VapourDiffusion/CreateVapourDiffusionFEBEX.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
+ *
+ * Created on March 5, 2021, 4:51 PM
+ */
+
+#include "CreateVapourDiffusionFEBEX.h"
+
+#include "BaseLib/ConfigTree.h"
+#include "MaterialLib/MPL/Property.h"
+#include "VapourDiffusionFEBEX.h"
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<Property> createVapourDiffusionFEBEX(
+    BaseLib::ConfigTree const& config)
+{
+    //! \ogs_file_param{properties__property__type}
+    config.checkConfigParameter("type", "VapourDiffusionFEBEX");
+    DBUG("Create VapourDiffusionFEBEX medium property");
+
+    // Second access for storage.
+    //! \ogs_file_param{properties__property__name}
+    auto property_name = config.peekConfigParameter<std::string>("name");
+
+    auto const tortuosity =
+        //! \ogs_file_param{properties__property__VapourDiffusionFEBEX__tortuosity}
+        config.getConfigParameter<double>("tortuosity");
+
+    return std::make_unique<VapourDiffusionFEBEX>(std::move(property_name),
+                                                  tortuosity);
+}
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/VapourDiffusion/CreateVapourDiffusionFEBEX.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
+ *
+ * Created on March 5, 2021, 4:51 PM
+ */
+
+#pragma once
+
+#include <memory>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace MaterialPropertyLib
+{
+class Property;
+std::unique_ptr<Property> createVapourDiffusionFEBEX(
+    BaseLib::ConfigTree const& config);
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionFEBEX.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
+ *
+ * Created on March 5, 2021, 3:49 PM
+ */
+
+#include "VapourDiffusionFEBEX.h"
+
+#include <algorithm>
+#include <cmath>
+
+#include "MaterialLib/MPL/Medium.h"
+#include "MaterialLib/MPL/VariableType.h"
+#include "MaterialLib/PhysicalConstant.h"
+
+namespace MaterialPropertyLib
+{
+PropertyDataType VapourDiffusionFEBEX::value(
+    const VariableArray& variable_array,
+    const ParameterLib::SpatialPosition& /*pos*/, const double /*t*/,
+    const double /*dt*/) const
+{
+    const double S_L = std::clamp(
+        std::get<double>(
+            variable_array[static_cast<int>(Variable::liquid_saturation)]),
+        0.0, 1.0);
+
+    const double T = std::get<double>(
+        variable_array[static_cast<int>(Variable::temperature)]);
+    const double phi =
+        std::get<double>(variable_array[static_cast<int>(Variable::porosity)]);
+
+    const double D_vr = tortuosity_ * phi * (1 - S_L);
+
+    return 2.16e-5 *
+           std::pow(T / MaterialLib::PhysicalConstant::CelsiusZeroInKelvin,
+                    1.8) *
+           D_vr;
+}
+
+PropertyDataType VapourDiffusionFEBEX::dValue(
+    VariableArray const& variable_array, Variable const variable,
+    ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
+    double const /*dt*/) const
+{
+    const double S_L = std::clamp(
+        std::get<double>(
+            variable_array[static_cast<int>(Variable::liquid_saturation)]),
+        0.0, 1.0);
+
+    const double T = std::get<double>(
+        variable_array[static_cast<int>(Variable::temperature)]);
+    const double phi =
+        std::get<double>(variable_array[static_cast<int>(Variable::porosity)]);
+
+    if (variable == Variable::temperature)
+    {
+        const double D_vr = tortuosity_ * phi * (1 - S_L);
+
+        return 1.8 * 2.16e-5 *
+               std::pow(T / MaterialLib::PhysicalConstant::CelsiusZeroInKelvin,
+                        0.8) *
+               D_vr / MaterialLib::PhysicalConstant::CelsiusZeroInKelvin;
+    }
+    if (variable == Variable::liquid_saturation)
+    {
+        return -2.16e-5 *
+               std::pow(T / MaterialLib::PhysicalConstant::CelsiusZeroInKelvin,
+                        1.8) *
+               tortuosity_ * phi;
+    }
+
+    OGS_FATAL(
+        "VapourDiffusionFEBEX::dValue is implemented for "
+        "derivatives with respect to temperature or saturation only.");
+}
+
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionFEBEX.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
+ *
+ * Created on March 5, 2021, 3:49 PM
+ */
+
+#pragma once
+
+#include "MaterialLib/MPL/Property.h"
+
+namespace MaterialPropertyLib
+{
+class Phase;
+
+/**
+ * \brief FEBEX type Vapour diffusion
+ *
+ *  The model was presented in \cite Rutquist2007TaskA1.
+ *
+ *  The vapour diffusion can be described by
+ *   \cite moldrup1997modeling, \cite moldrup2000predicting,
+ *  \cite chau2005simulation
+ *  \f[
+ *     D_v=2.16\cdot 10^{-5} \left(\frac{T}{273.15}\right)^{1.8}
+ *     D_{vr},
+ *  \f]
+ *  where \f$D_{vr}\f$ is the the relative diffusion coefficient,
+ *  and \f$T\f$ is the temperature.
+ *
+ *  In the FEBEX type, \f$D_{vr}\f$ takes the form of \cite Rutquist2007TaskA1
+ *   \f[
+ *     D_{vr}=\tau \phi (1 - S ),
+ *   \f]
+ *    with \f$\phi\f$, the porosity, \f$S\f$, the water saturation,
+ *    and  \f$\tau\f$ the tortuosity.
+ *
+ */
+class VapourDiffusionFEBEX final : public Property
+{
+public:
+    VapourDiffusionFEBEX(std::string name, double const tortuosity)
+        : tortuosity_(tortuosity)
+    {
+        name_ = std::move(name);
+    }
+
+    void checkScale() const override
+    {
+        if (!std::holds_alternative<Phase*>(scale_))
+        {
+            OGS_FATAL(
+                "The property 'VapourDiffusionFEBEX' is "
+                "implemented on the 'phase' scale only.");
+        }
+    }
+
+    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:
+    double const tortuosity_;
+};
+
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/PropertyType.h

@@ -94,6 +94,7 @@ enum PropertyType : int
     /// used to compute the hydrodynamic dispersion tensor.
     transversal_dispersivity,
     vapor_pressure,
+    vapour_diffusion,
     viscosity,
     volume_fraction,
     number_of_properties
@@ -314,6 +315,10 @@ inline PropertyType convertStringToProperty(std::string const& inString)
     {
         return PropertyType::vapor_pressure;
     }
+    if (boost::iequals(inString, "vapour_diffusion"))
+    {
+        return PropertyType::vapour_diffusion;
+    }
     if (boost::iequals(inString, "viscosity"))
     {
         return PropertyType::viscosity;
@@ -384,6 +389,7 @@ static const std::array<std::string, PropertyType::number_of_properties>
                              "transport_porosity",
                              "transversal_dispersivity",
                              "vapor_pressure",
+                             "vapour_diffusion",
                              "viscosity",
                              "volume_fraction"}};
 

Tests/MaterialLib/TestMPLVapourDiffusionFEBEX.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
+ *
+ * Created on March 4, 2021, 5:23 PM
+ */
+
+#include <gtest/gtest.h>
+
+#include <array>
+#include <cmath>
+#include <limits>
+
+#include "MaterialLib/MPL/Properties/VapourDiffusion/CreateVapourDiffusionFEBEX.h"
+#include "TestMPL.h"
+
+TEST(MaterialPropertyLib, VapourDiffusionFEBEX)
+{
+    char const xml[] =
+        "<property>"
+        "   <name>vapour_diffusion</name>"
+        "   <type>VapourDiffusionFEBEX</type>"
+        "   <tortuosity>0.8</tortuosity>"
+        "</property>";
+
+    std::unique_ptr<MaterialPropertyLib::Property> const property_ptr =
+        Tests::createTestProperty(
+            xml, MaterialPropertyLib::createVapourDiffusionFEBEX);
+    MaterialPropertyLib::Property const& property = *property_ptr;
+
+    double const T = 290.0;
+    double const S = 0.5;
+    double const phi = 0.15;
+
+    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();
+    variable_array[static_cast<int>(
+        MaterialPropertyLib::Variable::temperature)] = T;
+    variable_array[static_cast<int>(
+        MaterialPropertyLib::Variable::liquid_saturation)] = S;
+    variable_array[static_cast<int>(MaterialPropertyLib::Variable::porosity)] =
+        phi;
+
+    // The derivative of the water vapour with respect of temperature
+    {
+        std::array const Ts = {273.0, 293.0, 393.0, 420.0, 500.0};
+        std::array const D_v_expected = {1.294719e-06, 1.470431e-06,
+                                         2.494534e-06, 2.811458e-06,
+                                         3.847945e-06};
+
+        for (std::size_t i = 0; i < Ts.size(); ++i)
+        {
+            variable_array[static_cast<int>(
+                MaterialPropertyLib::Variable::temperature)] = Ts[i];
+
+            double const D_v =
+                property.template value<double>(variable_array, pos, t, dt);
+
+            ASSERT_LE(std::fabs(D_v_expected[i] - D_v), 1e-10)
+                << "for expected water vapour diffusion " << D_v_expected[i]
+                << " and for computed water vapour diffusion " << D_v;
+
+            double const dT = 1.0e-4;
+            variable_array[static_cast<int>(
+                MaterialPropertyLib::Variable::temperature)] = Ts[i] - dT;
+            double const D_v0 =
+                property.template value<double>(variable_array, pos, t, dt);
+
+            variable_array[static_cast<int>(
+                MaterialPropertyLib::Variable::temperature)] = Ts[i] + dT;
+            double const D_v1 =
+                property.template value<double>(variable_array, pos, t, dt);
+
+            double const approximated_dDv_dT = 0.5 * (D_v1 - D_v0) / dT;
+
+            double const analytic_dDv_dT = property.template dValue<double>(
+                variable_array, MaterialPropertyLib::Variable::temperature, pos,
+                t, dt);
+
+            ASSERT_LE(std::fabs(approximated_dDv_dT - analytic_dDv_dT), 1e-7)
+                << "for expected derivative of water vapour diffusion with "
+                   "respect to temperature "
+                << approximated_dDv_dT
+                << " and for computed derivative of water vapour diffusion "
+                   "with respect to temperature."
+                << analytic_dDv_dT;
+        }
+    }
+    // The derivative of the water vapour with respect of saturation
+    {
+        std::array const S = {-1.0, 0.0, 0.2,  0.33, 0.45,
+                              0.52, 0.6, 0.85, 1.0,  1.1};
+        std::array const D_v_expected = {
+            2.886883e-06, 2.886883e-06, 2.309507e-06, 1.934212e-06,
+            1.587786e-06, 1.385704e-06, 1.154753e-06, 4.330325e-07,
+            0.000000e+00, 0.000000e+00};
+        for (std::size_t i = 0; i < S.size(); ++i)
+        {
+            variable_array[static_cast<int>(
+                MaterialPropertyLib::Variable::temperature)] = T;
+
+            double const S_L_i = S[i];
+            variable_array[static_cast<int>(MPL::Variable::liquid_saturation)] =
+                S_L_i;
+            double const D_v =
+                property.template value<double>(variable_array, pos, t, dt);
+
+            ASSERT_LE(std::fabs(D_v_expected[i] - D_v), 1e-10)
+                << "for expected water vapour diffusion " << D_v_expected[i]
+                << " and for computed water vapour diffusion " << D_v;
+
+            double const analytic_dDv_dS = property.template dValue<double>(
+                variable_array, MPL::Variable::liquid_saturation, pos, t, dt);
+            double const dS_L = 1e-8;
+            double S_L_a = S_L_i - dS_L;
+            double S_L_b = S_L_i + dS_L;
+            double factor = 0.5;
+
+            if (S_L_i <= 0.0)
+            {
+                S_L_a = S_L_i;
+                S_L_b = dS_L;
+                factor = 1.0;
+            }
+            else if (S_L_i >= 1.0)
+            {
+                S_L_a = 1.0 - dS_L;
+                S_L_b = S_L_i;
+                factor = 1.0;
+            }
+
+            variable_array[static_cast<int>(MPL::Variable::liquid_saturation)] =
+                S_L_a;
+            double const D_v_a =
+                property.template value<double>(variable_array, pos, t, dt);
+            variable_array[static_cast<int>(MPL::Variable::liquid_saturation)] =
+                S_L_b;
+            double const D_v_b =
+                property.template value<double>(variable_array, pos, t, dt);
+
+            double const approximated_dDv_dS = factor * (D_v_b - D_v_a) / dS_L;
+
+            ASSERT_LE(std::fabs(approximated_dDv_dS - analytic_dDv_dS) /
+                          analytic_dDv_dS,
+                      1e-10)
+                << "for expected derivative of water vapour diffusion with "
+                   "respect to saturation "
+                << approximated_dDv_dS
+                << " and for computed derivative of water vapour diffusion "
+                   "with respect to saturation."
+                << analytic_dDv_dS;
+        }
+    }
+}