From a071eaf01b02b2c0e287bf4a9cf087315ba06cd4 Mon Sep 17 00:00:00 2001 From: Wenqing Wang Date: Fri, 16 Jun 2023 15:22:54 +0200 Subject: [PATCH 1/6] [MPL] Removed n(1-S_L) from VapourDiffusionPMQ --- .../VapourDiffusion/VapourDiffusionPMQ.cpp | 8 +++--- .../VapourDiffusion/VapourDiffusionPMQ.h | 27 +++++++++++++------ .../MaterialLib/TestMPLVapourDiffusionPMQ.cpp | 13 ++++----- 3 files changed, 27 insertions(+), 21 deletions(-) diff --git a/MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionPMQ.cpp b/MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionPMQ.cpp index b07f3bd4082..67bcfd6ec2f 100644 --- a/MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionPMQ.cpp +++ b/MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionPMQ.cpp @@ -28,10 +28,9 @@ 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 * std::pow(T / MaterialLib::PhysicalConstant::CelsiusZeroInKelvin, @@ -47,12 +46,11 @@ 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 * std::pow(T / MaterialLib::PhysicalConstant::CelsiusZeroInKelvin, @@ -65,7 +63,7 @@ PropertyDataType VapourDiffusionPMQ::dValue( return -2.16e-5 * std::pow(T / MaterialLib::PhysicalConstant::CelsiusZeroInKelvin, 1.8) * - 2.0 * tortuosity_ * phi * S_L; + tortuosity_; } OGS_FATAL( diff --git a/MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionPMQ.h b/MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionPMQ.h index 7e593d872d7..6dea2205f1a 100644 --- a/MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionPMQ.h +++ b/MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionPMQ.h @@ -20,11 +20,8 @@ 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_{vr}, @@ -32,15 +29,29 @@ class Phase; * where \f$D_{vr}\f$ is the the relative diffusion coefficient, * and \f$T\f$ is the temperature. * - * In the PMQ type, \f$D_{vr}\f$ takes the form of \cite chau2005simulation + * 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. + * + * 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: diff --git a/Tests/MaterialLib/TestMPLVapourDiffusionPMQ.cpp b/Tests/MaterialLib/TestMPLVapourDiffusionPMQ.cpp index 1770f4a3a8b..4982f8a4257 100644 --- a/Tests/MaterialLib/TestMPLVapourDiffusionPMQ.cpp +++ b/Tests/MaterialLib/TestMPLVapourDiffusionPMQ.cpp @@ -33,7 +33,6 @@ TEST(MaterialPropertyLib, VapourDiffusionPMQ) const double T = 290.0; const double S = 0.5; - const double phi = 0.15; MaterialPropertyLib::VariableArray variable_array; ParameterLib::SpatialPosition const pos; @@ -41,14 +40,13 @@ TEST(MaterialPropertyLib, VapourDiffusionPMQ) double const dt = std::numeric_limits::quiet_NaN(); variable_array.temperature = T; variable_array.liquid_saturation = S; - variable_array.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 = {5.340717e-07, 6.065526e-07, - 1.028995e-06, 1.159726e-06, - 1.587277e-06}; + std::array const D_v_expected = { + 7.1209560000000006e-06, 8.087368e-06, 1.3719933333333336e-05, + 1.5463013333333333e-05, 2.1163693333333336e-05}; for (std::size_t i = 0; i < Ts.size(); ++i) { @@ -90,9 +88,8 @@ TEST(MaterialPropertyLib, VapourDiffusionPMQ) 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.381679e-06, 2.381679e-06, 1.524274e-06, 1.069136e-06, - 7.204578e-07, 5.487388e-07, 3.810686e-07, 5.358777e-08, - 0.000000e+00, 0.000000e+00}; + 1.58779e-05, 1.58779e-05, 1.27023e-05, 1.06382e-05, 8.73282e-06, + 7.62137e-06, 6.35114e-06, 2.38168e-06, 0.0, 0.0}; for (std::size_t i = 0; i < S.size(); ++i) { variable_array.temperature = T; From 5b2b991acd1b3632b41a4ad8add28f4e8eb96f9e Mon Sep 17 00:00:00 2001 From: Wenqing Wang Date: Fri, 16 Jun 2023 16:26:55 +0200 Subject: [PATCH 2/6] [MPL] Removed n(1-S_L) from VapourDiffusionFEBEX --- .../VapourDiffusion/VapourDiffusionFEBEX.cpp | 19 +------ .../VapourDiffusion/VapourDiffusionFEBEX.h | 10 +++- .../TestMPLVapourDiffusionFEBEX.cpp | 55 ++++--------------- 3 files changed, 20 insertions(+), 64 deletions(-) diff --git a/MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionFEBEX.cpp b/MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionFEBEX.cpp index 94f66033439..b9488d5ba74 100644 --- a/MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionFEBEX.cpp +++ b/MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionFEBEX.cpp @@ -25,17 +25,12 @@ 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 * std::pow(T / MaterialLib::PhysicalConstant::CelsiusZeroInKelvin, 1.8) * - D_vr; + tortuosity_; } PropertyDataType VapourDiffusionFEBEX::dValue( @@ -43,26 +38,18 @@ PropertyDataType VapourDiffusionFEBEX::dValue( 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 * std::pow(T / MaterialLib::PhysicalConstant::CelsiusZeroInKelvin, 0.8) * - D_vr / MaterialLib::PhysicalConstant::CelsiusZeroInKelvin; + 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( diff --git a/MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionFEBEX.h b/MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionFEBEX.h index 7e015c2d49b..b52f6fd958a 100644 --- a/MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionFEBEX.h +++ b/MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionFEBEX.h @@ -34,11 +34,15 @@ class Phase; * * 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 { diff --git a/Tests/MaterialLib/TestMPLVapourDiffusionFEBEX.cpp b/Tests/MaterialLib/TestMPLVapourDiffusionFEBEX.cpp index 8bdabb4691b..011bd7536a0 100644 --- a/Tests/MaterialLib/TestMPLVapourDiffusionFEBEX.cpp +++ b/Tests/MaterialLib/TestMPLVapourDiffusionFEBEX.cpp @@ -33,23 +33,18 @@ TEST(MaterialPropertyLib, VapourDiffusionFEBEX) 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::quiet_NaN(); double const dt = std::numeric_limits::quiet_NaN(); variable_array.temperature = T; - variable_array.liquid_saturation = S; - variable_array.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}; + std::array const D_v_expected = {1.72629e-05, 1.96057e-05, 3.32605e-05, + 3.74861e-05, 5.13059e-05, 1.92459e-05}; for (std::size_t i = 0; i < Ts.size(); ++i) { @@ -86,62 +81,32 @@ TEST(MaterialPropertyLib, VapourDiffusionFEBEX) << analytic_dDv_dT; } } - // The derivative of the water vapour with respect of saturation + // The derivative of the water vapour with respect of saturation, which is + // zero. { 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}; + double const D_v_expected = 1.92459e-05; for (std::size_t i = 0; i < S.size(); ++i) { variable_array.temperature = T; - double const S_L_i = S[i]; - variable_array.liquid_saturation = S_L_i; double const D_v = property.template value(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] + ASSERT_LE(std::fabs(D_v_expected - D_v), 1e-10) + << "for expected water vapour diffusion " << D_v_expected << " and for computed water vapour diffusion " << D_v; double const analytic_dDv_dS = property.template dValue( 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.liquid_saturation = S_L_a; - double const D_v_a = - property.template value(variable_array, pos, t, dt); - variable_array.liquid_saturation = S_L_b; - double const D_v_b = - property.template value(variable_array, pos, t, dt); - double const approximated_dDv_dS = factor * (D_v_b - D_v_a) / dS_L; + double const expected_dDv_dS = 0.0; - ASSERT_LE(std::fabs(approximated_dDv_dS - analytic_dDv_dS) / - analytic_dDv_dS, - 1e-10) + ASSERT_LE(std::fabs(expected_dDv_dS - analytic_dDv_dS), 1e-10) << "for expected derivative of water vapour diffusion with " "respect to saturation " - << approximated_dDv_dS + << expected_dDv_dS << " and for computed derivative of water vapour diffusion " "with respect to saturation." << analytic_dDv_dS; From fcfc828dc9f75ef5f213035b64b8c5fe8c838215 Mon Sep 17 00:00:00 2001 From: Wenqing Wang Date: Fri, 16 Jun 2023 16:49:41 +0200 Subject: [PATCH 3/6] [ProcessLib] Update the usage of the vapour diffusion model --- .../ThermoRichardsFlow/ThermoRichardsFlowFEM-impl.h | 2 ++ .../ConstitutiveCommon/TRMVaporDiffusion.cpp | 13 +++++++------ 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowFEM-impl.h b/ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowFEM-impl.h index c3dd6ccb814..f3805644612 100644 --- a/ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowFEM-impl.h +++ b/ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowFEM-impl.h @@ -534,6 +534,7 @@ void ThermoRichardsFlowLocalAssembler:: variables.porosity = phi; double const D_v = + phi * (1.0 - S_L) * gas_phase->property(MPL::PropertyType::diffusion) .template value(variables, x_position, t, dt); @@ -990,6 +991,7 @@ void ThermoRichardsFlowLocalAssembler::assemble( variables.porosity = phi; double const D_v = + phi * (1.0 - S_L) * gas_phase->property(MPL::PropertyType::diffusion) .template value(variables, x_position, t, dt); diff --git a/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/TRMVaporDiffusion.cpp b/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/TRMVaporDiffusion.cpp index 2435e1ea94e..a7c17a5430d 100644 --- a/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/TRMVaporDiffusion.cpp +++ b/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/TRMVaporDiffusion.cpp @@ -74,7 +74,10 @@ void TRMVaporDiffusionModel::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(variables, x_t.x, x_t.t, x_t.dt); @@ -89,21 +92,19 @@ void TRMVaporDiffusionModel::eval( MaterialPropertyLib::PropertyType::specific_heat_capacity) .template value(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) From 27eb293555b85406ccd507795dc2c054f375b046 Mon Sep 17 00:00:00 2001 From: Wenqing Wang Date: Mon, 4 Sep 2023 14:41:23 +0200 Subject: [PATCH 4/6] [MPL] Add two optional input tags for the base diffusion coefficient and exponent, respectively, of vapour diffusion models --- .../CreateVapourDiffusionFEBEX.cpp | 14 +++++++++++-- .../CreateVapourDiffusionPMQ.cpp | 12 ++++++++++- .../VapourDiffusion/VapourDiffusionFEBEX.cpp | 8 +++---- .../VapourDiffusion/VapourDiffusionFEBEX.h | 20 +++++++++++++----- .../VapourDiffusion/VapourDiffusionPMQ.cpp | 12 +++++------ .../VapourDiffusion/VapourDiffusionPMQ.h | 21 +++++++++++++++---- 6 files changed, 65 insertions(+), 22 deletions(-) diff --git a/MaterialLib/MPL/Properties/VapourDiffusion/CreateVapourDiffusionFEBEX.cpp b/MaterialLib/MPL/Properties/VapourDiffusion/CreateVapourDiffusionFEBEX.cpp index 4f65c5ac862..7d9d899bdb1 100644 --- a/MaterialLib/MPL/Properties/VapourDiffusion/CreateVapourDiffusionFEBEX.cpp +++ b/MaterialLib/MPL/Properties/VapourDiffusion/CreateVapourDiffusionFEBEX.cpp @@ -32,7 +32,17 @@ std::unique_ptr createVapourDiffusionFEBEX( //! \ogs_file_param{properties__property__VapourDiffusionFEBEX__tortuosity} config.getConfigParameter("tortuosity"); - return std::make_unique(std::move(property_name), - tortuosity); + double const base_diffusion_coefficient = + //! \ogs_file_param{properties__property__VapourDiffusionFEBEX__base_diffusion_coefficient} + config.getConfigParameter("base_diffusion_coefficient", + 2.16e-5); + + double const exponent = + //! \ogs_file_param{properties__property__VapourDiffusionFEBEX__exponent} + config.getConfigParameter("exponent", 1.8); + + return std::make_unique( + std::move(property_name), tortuosity, base_diffusion_coefficient, + exponent); } } // namespace MaterialPropertyLib diff --git a/MaterialLib/MPL/Properties/VapourDiffusion/CreateVapourDiffusionPMQ.cpp b/MaterialLib/MPL/Properties/VapourDiffusion/CreateVapourDiffusionPMQ.cpp index a303fcc828d..fe9ff80e7e5 100644 --- a/MaterialLib/MPL/Properties/VapourDiffusion/CreateVapourDiffusionPMQ.cpp +++ b/MaterialLib/MPL/Properties/VapourDiffusion/CreateVapourDiffusionPMQ.cpp @@ -27,7 +27,17 @@ std::unique_ptr createVapourDiffusionPMQ( //! \ogs_file_param{properties__property__name} auto property_name = config.peekConfigParameter("name"); + double const base_diffusion_coefficient = + //! \ogs_file_param{properties__property__VapourDiffusionPMQ__base_diffusion_coefficient} + config.getConfigParameter("base_diffusion_coefficient", + 2.16e-5); + + double const exponent = + //! \ogs_file_param{properties__property__VapourDiffusionPMQ__exponent} + config.getConfigParameter("exponent", 1.8); + //! \ogs_file_param_special{properties__property__VapourDiffusionPMQ} - return std::make_unique(std::move(property_name)); + return std::make_unique( + std::move(property_name), base_diffusion_coefficient, exponent); } } // namespace MaterialPropertyLib diff --git a/MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionFEBEX.cpp b/MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionFEBEX.cpp index b9488d5ba74..8f46f8be6ce 100644 --- a/MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionFEBEX.cpp +++ b/MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionFEBEX.cpp @@ -27,9 +27,9 @@ PropertyDataType VapourDiffusionFEBEX::value( { const double T = variable_array.temperature; - return 2.16e-5 * + return base_diffusion_coefficient_ * std::pow(T / MaterialLib::PhysicalConstant::CelsiusZeroInKelvin, - 1.8) * + exponent_) * tortuosity_; } @@ -42,9 +42,9 @@ PropertyDataType VapourDiffusionFEBEX::dValue( if (variable == Variable::temperature) { - return 1.8 * 2.16e-5 * + return exponent_ * base_diffusion_coefficient_ * std::pow(T / MaterialLib::PhysicalConstant::CelsiusZeroInKelvin, - 0.8) * + exponent_ - 1.0) * tortuosity_ / MaterialLib::PhysicalConstant::CelsiusZeroInKelvin; } if (variable == Variable::liquid_saturation) diff --git a/MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionFEBEX.h b/MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionFEBEX.h index b52f6fd958a..0711b1b11ba 100644 --- a/MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionFEBEX.h +++ b/MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionFEBEX.h @@ -26,11 +26,15 @@ 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[ @@ -47,8 +51,12 @@ class Phase; 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); } @@ -75,6 +83,8 @@ class VapourDiffusionFEBEX final : public Property private: double const tortuosity_; + double const base_diffusion_coefficient_; + double const exponent_; }; } // namespace MaterialPropertyLib diff --git a/MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionPMQ.cpp b/MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionPMQ.cpp index 67bcfd6ec2f..ea6ff097b10 100644 --- a/MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionPMQ.cpp +++ b/MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionPMQ.cpp @@ -32,9 +32,9 @@ PropertyDataType VapourDiffusionPMQ::value( const double S_v = 1 - S_L; 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; } @@ -52,17 +52,17 @@ PropertyDataType VapourDiffusionPMQ::dValue( const double S_v = 1 - S_L; 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) * + exponent_) * tortuosity_; } diff --git a/MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionPMQ.h b/MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionPMQ.h index 6dea2205f1a..692dd1edef5 100644 --- a/MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionPMQ.h +++ b/MaterialLib/MPL/Properties/VapourDiffusion/VapourDiffusionPMQ.h @@ -23,11 +23,15 @@ class Phase; * The vapour diffusion can be described by * \cite moldrup1997modeling, \cite moldrup2000predicting, * \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 @@ -55,7 +59,14 @@ class Phase; 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 { @@ -79,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 From 17ba18cfec464988b5b28cc661e2e7b5b06d81d6 Mon Sep 17 00:00:00 2001 From: Wenqing Wang Date: Mon, 4 Sep 2023 14:42:50 +0200 Subject: [PATCH 5/6] [UnitTest/MPL] Test the newly added optional tag for VapourDiffusionFEBEX --- Tests/MaterialLib/TestMPLVapourDiffusionFEBEX.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Tests/MaterialLib/TestMPLVapourDiffusionFEBEX.cpp b/Tests/MaterialLib/TestMPLVapourDiffusionFEBEX.cpp index 011bd7536a0..0c1aa4b1b51 100644 --- a/Tests/MaterialLib/TestMPLVapourDiffusionFEBEX.cpp +++ b/Tests/MaterialLib/TestMPLVapourDiffusionFEBEX.cpp @@ -25,6 +25,8 @@ TEST(MaterialPropertyLib, VapourDiffusionFEBEX) " vapour_diffusion" " VapourDiffusionFEBEX" " 0.8" + " 2.16e-5" + " 1.8" ""; std::unique_ptr const property_ptr = From 8b85947471608eefc67e85891ef06feefeeebafa Mon Sep 17 00:00:00 2001 From: Wenqing Wang Date: Mon, 4 Sep 2023 14:51:46 +0200 Subject: [PATCH 6/6] [doc] Parameter documentation for the newly added input tag for base diffusion coefficient --- .../VapourDiffusionFEBEX/t_base_diffusion_coefficient.md | 3 +++ .../properties/property/VapourDiffusionFEBEX/t_exponent.md | 1 + .../VapourDiffusionPMQ/t_base_diffusion_coefficient.md | 3 +++ .../properties/property/VapourDiffusionPMQ/t_exponent.md | 1 + 4 files changed, 8 insertions(+) create mode 100644 Documentation/ProjectFile/properties/property/VapourDiffusionFEBEX/t_base_diffusion_coefficient.md create mode 100644 Documentation/ProjectFile/properties/property/VapourDiffusionFEBEX/t_exponent.md create mode 100644 Documentation/ProjectFile/properties/property/VapourDiffusionPMQ/t_base_diffusion_coefficient.md create mode 100644 Documentation/ProjectFile/properties/property/VapourDiffusionPMQ/t_exponent.md diff --git a/Documentation/ProjectFile/properties/property/VapourDiffusionFEBEX/t_base_diffusion_coefficient.md b/Documentation/ProjectFile/properties/property/VapourDiffusionFEBEX/t_base_diffusion_coefficient.md new file mode 100644 index 00000000000..90d05acc767 --- /dev/null +++ b/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. diff --git a/Documentation/ProjectFile/properties/property/VapourDiffusionFEBEX/t_exponent.md b/Documentation/ProjectFile/properties/property/VapourDiffusionFEBEX/t_exponent.md new file mode 100644 index 00000000000..4299bfcf5a0 --- /dev/null +++ b/Documentation/ProjectFile/properties/property/VapourDiffusionFEBEX/t_exponent.md @@ -0,0 +1 @@ +An optional input to the exponent, which has a default value 1.8. diff --git a/Documentation/ProjectFile/properties/property/VapourDiffusionPMQ/t_base_diffusion_coefficient.md b/Documentation/ProjectFile/properties/property/VapourDiffusionPMQ/t_base_diffusion_coefficient.md new file mode 100644 index 00000000000..90d05acc767 --- /dev/null +++ b/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. diff --git a/Documentation/ProjectFile/properties/property/VapourDiffusionPMQ/t_exponent.md b/Documentation/ProjectFile/properties/property/VapourDiffusionPMQ/t_exponent.md new file mode 100644 index 00000000000..4299bfcf5a0 --- /dev/null +++ b/Documentation/ProjectFile/properties/property/VapourDiffusionPMQ/t_exponent.md @@ -0,0 +1 @@ +An optional input to the exponent, which has a default value 1.8.