Merge branch 'impr_vapor_diffusion' into 'master'

[MPL] removed phi (1-S_L) from the vapour diffusion model

See merge request ogs/ogs!4648

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

@@ -0,0 +1,3 @@
+An optional input to the base diffusion coefficient, which has a default value
+ 2.16e-5 \f${\text m}^2 \text{Pa}/(\text{s}\text{K}^{n})\f$ with \f$n\f$, the
+ exponent.

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

@@ -0,0 +1 @@
+An optional input to the exponent, which has a default value 1.8.

Documentation/ProjectFile/properties/property/VapourDiffusionPMQ/t_base_diffusion_coefficient.md

@@ -0,0 +1,3 @@
+An optional input to the base diffusion coefficient, which has a default value
+ 2.16e-5 \f${\text m}^2 \text{Pa}/(\text{s}\text{K}^{n})\f$ with \f$n\f$, the
+ exponent.

Documentation/ProjectFile/properties/property/VapourDiffusionPMQ/t_exponent.md

@@ -0,0 +1 @@
+An optional input to the exponent, which has a default value 1.8.

MaterialLib/MPL/Properties/VapourDiffusion/CreateVapourDiffusionFEBEX.cpp

@@ -32,7 +32,17 @@ std::unique_ptr<Property> createVapourDiffusionFEBEX(
         //! \ogs_file_param{properties__property__VapourDiffusionFEBEX__tortuosity}
         config.getConfigParameter<double>("tortuosity");
 
-    return std::make_unique<VapourDiffusionFEBEX>(std::move(property_name),
-                                                  tortuosity);
+    double const base_diffusion_coefficient =
+        //! \ogs_file_param{properties__property__VapourDiffusionFEBEX__base_diffusion_coefficient}
+        config.getConfigParameter<double>("base_diffusion_coefficient",
+                                          2.16e-5);
+
+    double const exponent =
+        //! \ogs_file_param{properties__property__VapourDiffusionFEBEX__exponent}
+        config.getConfigParameter<double>("exponent", 1.8);
+
+    return std::make_unique<VapourDiffusionFEBEX>(
+        std::move(property_name), tortuosity, base_diffusion_coefficient,
+        exponent);
 }
 }  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/VapourDiffusion/CreateVapourDiffusionPMQ.cpp

@@ -27,7 +27,17 @@ std::unique_ptr<Property> createVapourDiffusionPMQ(
     //! \ogs_file_param{properties__property__name}
     auto property_name = config.peekConfigParameter<std::string>("name");
 
+    double const base_diffusion_coefficient =
+        //! \ogs_file_param{properties__property__VapourDiffusionPMQ__base_diffusion_coefficient}
+        config.getConfigParameter<double>("base_diffusion_coefficient",
+                                          2.16e-5);
+
+    double const exponent =
+        //! \ogs_file_param{properties__property__VapourDiffusionPMQ__exponent}
+        config.getConfigParameter<double>("exponent", 1.8);
+
     //! \ogs_file_param_special{properties__property__VapourDiffusionPMQ}
-    return std::make_unique<VapourDiffusionPMQ>(std::move(property_name));
+    return std::make_unique<VapourDiffusionPMQ>(
+        std::move(property_name), base_diffusion_coefficient, exponent);
 }
 }  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionFEBEX.cpp

@@ -25,44 +25,31 @@ PropertyDataType VapourDiffusionFEBEX::value(
     const ParameterLib::SpatialPosition& /*pos*/, const double /*t*/,
     const double /*dt*/) const
 {
-    const double S_L = std::clamp(variable_array.liquid_saturation, 0.0, 1.0);
-
     const double T = variable_array.temperature;
-    const double phi = variable_array.porosity;
-
-    const double D_vr = tortuosity_ * phi * (1 - S_L);
 
-    return 2.16e-5 *
+    return base_diffusion_coefficient_ *
            std::pow(T / MaterialLib::PhysicalConstant::CelsiusZeroInKelvin,
-                    1.8) *
-           D_vr;
+                    exponent_) *
+           tortuosity_;
 }
 
 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(variable_array.liquid_saturation, 0.0, 1.0);
-
     const double T = variable_array.temperature;
-    const double phi = variable_array.porosity;
 
     if (variable == Variable::temperature)
     {
-        const double D_vr = tortuosity_ * phi * (1 - S_L);
-
-        return 1.8 * 2.16e-5 *
+        return exponent_ * base_diffusion_coefficient_ *
                std::pow(T / MaterialLib::PhysicalConstant::CelsiusZeroInKelvin,
-                        0.8) *
-               D_vr / MaterialLib::PhysicalConstant::CelsiusZeroInKelvin;
+                        exponent_ - 1.0) *
+               tortuosity_ / MaterialLib::PhysicalConstant::CelsiusZeroInKelvin;
     }
     if (variable == Variable::liquid_saturation)
     {
-        return -2.16e-5 *
-               std::pow(T / MaterialLib::PhysicalConstant::CelsiusZeroInKelvin,
-                        1.8) *
-               tortuosity_ * phi;
+        return 0.0;
     }
 
     OGS_FATAL(

MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionFEBEX.h

@@ -26,25 +26,37 @@ class Phase;
  *   \cite moldrup1997modeling, \cite moldrup2000predicting,
  *  \cite chau2005simulation
  *  \f[
- *     D_v=2.16\cdot 10^{-5} \left(\frac{T}{273.15}\right)^{1.8}
+ *     D_v=D_0 \left(\frac{T}{273.15}\right)^{n}
  *     D_{vr},
  *  \f]
- *  where \f$D_{vr}\f$ is the the relative diffusion coefficient,
- *  and \f$T\f$ is the temperature.
+ *  where \f$D_{0}\f$ is the base diffusion coefficient with default value
+ *  \f$2.16\cdot 10^{-5}\f$ \f${\text m}^2
+ * \text{Pa}/(\text{s}\text{K}^{n})\f$,
+ *  \f$n\f$ is the exponent with default value 1.8,
+ *  \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 ),
+ *     D_{vr}=\tau \phi (1-S_L),
  *   \f]
- *    with \f$\phi\f$, the porosity, \f$S\f$, the water saturation,
- *    and  \f$\tau\f$ the tortuosity.
+ *    with \f$\phi\f$, the porosity, \f$S_L\f$, the liquid saturation,
+ *    and  \f$\tau\f$, the tortuosity.
  *
+ *  Note: In order to maintain consistency with the implementation of the
+ *  computations of other vapor-related parameters, \f$ \phi (1 - S_L )\f$ is
+ *  removed from the implementation for this class and is multiplied back in
+ *  the local assembler.
  */
 class VapourDiffusionFEBEX final : public Property
 {
 public:
-    VapourDiffusionFEBEX(std::string name, double const tortuosity)
-        : tortuosity_(tortuosity)
+    VapourDiffusionFEBEX(std::string name, double const tortuosity,
+                         double const base_diffusion_coefficient,
+                         double const exponent)
+        : tortuosity_(tortuosity),
+          base_diffusion_coefficient_(base_diffusion_coefficient),
+          exponent_(exponent)
     {
         name_ = std::move(name);
     }
@@ -71,6 +83,8 @@ class VapourDiffusionFEBEX final : public Property
 
 private:
     double const tortuosity_;
+    double const base_diffusion_coefficient_;
+    double const exponent_;
 };
 
 }  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionPMQ.cpp

@@ -28,14 +28,13 @@ PropertyDataType VapourDiffusionPMQ::value(
     const double S_L = std::clamp(variable_array.liquid_saturation, 0.0, 1.0);
 
     const double T = variable_array.temperature;
-    const double phi = variable_array.porosity;
 
     const double S_v = 1 - S_L;
-    const double D_vr = tortuosity_ * phi * S_v * S_v;
+    const double D_vr = tortuosity_ * S_v;
 
-    return 2.16e-5 *
+    return base_diffusion_coefficient_ *
            std::pow(T / MaterialLib::PhysicalConstant::CelsiusZeroInKelvin,
-                    1.8) *
+                    exponent_) *
            D_vr;
 }
 
@@ -47,25 +46,24 @@ PropertyDataType VapourDiffusionPMQ::dValue(
     const double S_L = std::clamp(variable_array.liquid_saturation, 0.0, 1.0);
 
     const double T = variable_array.temperature;
-    const double phi = variable_array.porosity;
 
     if (variable == Variable::temperature)
     {
         const double S_v = 1 - S_L;
-        const double D_vr = tortuosity_ * phi * S_v * S_v;
+        const double D_vr = tortuosity_ * S_v;
 
-        return 1.8 * 2.16e-5 *
+        return exponent_ * base_diffusion_coefficient_ *
                std::pow(T / MaterialLib::PhysicalConstant::CelsiusZeroInKelvin,
-                        0.8) *
+                        exponent_ - 1.0) *
                D_vr / MaterialLib::PhysicalConstant::CelsiusZeroInKelvin;
     }
 
     if (variable == Variable::liquid_saturation)
     {
-        return -2.16e-5 *
+        return -base_diffusion_coefficient_ *
                std::pow(T / MaterialLib::PhysicalConstant::CelsiusZeroInKelvin,
-                        1.8) *
-               2.0 * tortuosity_ * phi * S_L;
+                        exponent_) *
+               tortuosity_;
     }
 
     OGS_FATAL(

MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionPMQ.h

@@ -20,31 +20,53 @@ class Phase;
 /**
  * \brief The Penman-Millington-Quirk (PMQ) Vapour diffusion model.
  *
- * The model was presented in \cite chau2005simulation
- *
  *  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_v=D_0 \left(\frac{T}{273.15}\right)^{n}
  * 	D_{vr},
  *  \f]
- *  where \f$D_{vr}\f$ is the the relative diffusion coefficient,
- *  and \f$T\f$ is the temperature.
+ *  where \f$D_{0}\f$ is the base diffusion coefficient with default value
+ *  \f$2.16\cdot 10^{-5}\f$ \f${\text m}^2
+ * \text{Pa}/(\text{s}\text{K}^{n})\f$,
+ *  \f$n\f$ is the exponent with default value 1.8,
+ *  \f$D_{vr}\f$ is the the relative
+ * diffusion coefficient, and \f$T\f$ is the temperature.
+ *
+ *  The Penman–Millington–Quirk (PMQ) model \cite moldrup1997modeling is given
+ *  as
+ *  \f[
+ *     D_{vr}=0.66 \phi \left(\frac{\kappa}{\phi}\right)^{\frac{12-m}{3}},
+ *  \f]
+ * where \f$\phi\f$ is the total porosity, \f$\kappa\f$ is the air filled
+ *  porosity, and \f$m\f$ is a fitting parameter. The air filled porosity is
+ *  defined as \f$ \kappa = \phi-\theta = \phi -S_L \phi \f$ with \f$\theta\f$
+ *  the liquid content, and \f$S_L\f$ the liquid saturation.
  *
- *  In the PMQ type, \f$D_{vr}\f$ takes the form of \cite chau2005simulation
+ * According to the study presented in \cite moldrup1997modeling, \f$m=6\f$
+ *  is the best fitting parameter for the sieved, repacked soils that the
+ *  authors tested. Therefore,  \f$m=6\f$ is used in the implementation, which
+ * gives
  *   \f[
- *     D_{vr}=0.66 \phi (1 - S_L )^2,
+ *     D_{vr}=0.66 \phi (1 - S_L )^2.
  *   \f]
- *    with \f$\phi\f$, the porosity, \f$S_L\f$, the water saturation, \f$0.66\f$
- *     the tortuosity constant.
  *
+ * Note: In order to maintain consistency with the implementation of the
+ *  computations of other vapor-related parameters, \f$ \phi (1 - S_L )\f$ is
+ *  removed from the implementation for this class and is multiplied back in
+ *  the local assembler.
  */
-
 class VapourDiffusionPMQ final : public Property
 {
 public:
-    explicit VapourDiffusionPMQ(std::string name) { name_ = std::move(name); }
+    explicit VapourDiffusionPMQ(std::string name,
+                                double const base_diffusion_coefficient,
+                                double const exponent)
+        : base_diffusion_coefficient_(base_diffusion_coefficient),
+          exponent_(exponent)
+    {
+        name_ = std::move(name);
+    }
 
     void checkScale() const override
     {
@@ -68,6 +90,8 @@ class VapourDiffusionPMQ final : public Property
 
 private:
     static double constexpr tortuosity_ = 0.66;
+    double const base_diffusion_coefficient_;
+    double const exponent_;
 };
 
 }  // namespace MaterialPropertyLib

ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowFEM-impl.h

@@ -534,6 +534,7 @@ void ThermoRichardsFlowLocalAssembler<ShapeFunction, GlobalDim>::
 
             variables.porosity = phi;
             double const D_v =
+                phi * (1.0 - S_L) *
                 gas_phase->property(MPL::PropertyType::diffusion)
                     .template value<double>(variables, x_position, t, dt);
 
@@ -990,6 +991,7 @@ void ThermoRichardsFlowLocalAssembler<ShapeFunction, GlobalDim>::assemble(
 
             variables.porosity = phi;
             double const D_v =
+                phi * (1.0 - S_L) *
                 gas_phase->property(MPL::PropertyType::diffusion)
                     .template value<double>(variables, x_position, t, dt);
 

ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/TRMVaporDiffusion.cpp

@@ -74,7 +74,10 @@ void TRMVaporDiffusionModel<DisplacementDim>::eval(
 
         double const phi = poro_data.phi;
         variables.porosity = phi;
+
+        double const S_g = 1.0 - S_L_data.S_L;
         double const D_v =
+            phi * S_g *
             gas_phase->property(MPL::PropertyType::diffusion)
                 .template value<double>(variables, x_t.x, x_t.t, x_t.dt);
 
@@ -89,21 +92,19 @@ void TRMVaporDiffusionModel<DisplacementDim>::eval(
                     MaterialPropertyLib::PropertyType::specific_heat_capacity)
                 .template value<double>(variables, x_t.x, x_t.t, x_t.dt);
 
-        out.M_TT_X_NTN +=
-            out.heat_capacity_vapor * rho_wv * (1 - S_L_data.S_L) * phi;
+        out.M_TT_X_NTN += out.heat_capacity_vapor * rho_wv * S_g * phi;
 
         out.storage_coefficient_by_water_vapor =
-            phi *
-            (rho_wv * dS_L_data.dS_L_dp_cap + (1 - S_L_data.S_L) * drho_wv_dp);
+            phi * (rho_wv * dS_L_data.dS_L_dp_cap + S_g * drho_wv_dp);
 
-        out.M_pT_X_NTN += phi * (1 - S_L_data.S_L) * drho_wv_dT;
+        out.M_pT_X_NTN += phi * S_g * drho_wv_dT;
 
         //
         // Latent heat term
         //
         if (gas_phase->hasProperty(MPL::PropertyType::specific_latent_heat))
         {
-            double const factor = phi * (1 - S_L_data.S_L) / rho_L_data.rho_LR;
+            double const factor = phi * S_g / rho_L_data.rho_LR;
             // The volumetric latent heat of vaporization of liquid water
             double const L0 =
                 gas_phase->property(MPL::PropertyType::specific_latent_heat)