Merge pull request #1613 from Yonghui56/pr_2phase_totalmassdensity_pl…

…_wip

Create two phase flow process with total mass density model

Applications/ApplicationsLib/ProjectData.cpp

@@ -35,16 +35,17 @@
 #include "ProcessLib/UncoupledProcessesTimeLoop.h"
 
 #include "ProcessLib/GroundwaterFlow/CreateGroundwaterFlowProcess.h"
+#include "ProcessLib/HT/CreateHTProcess.h"
 #include "ProcessLib/HeatConduction/CreateHeatConductionProcess.h"
 #include "ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.h"
 #include "ProcessLib/LIE/HydroMechanics/CreateHydroMechanicsProcess.h"
 #include "ProcessLib/LIE/SmallDeformation/CreateSmallDeformationProcess.h"
-#include "ProcessLib/RichardsFlow/CreateRichardsFlowProcess.h"
 #include "ProcessLib/LiquidFlow/CreateLiquidFlowProcess.h"
+#include "ProcessLib/RichardsFlow/CreateRichardsFlowProcess.h"
 #include "ProcessLib/SmallDeformation/CreateSmallDeformationProcess.h"
 #include "ProcessLib/TES/CreateTESProcess.h"
-#include "ProcessLib/HT/CreateHTProcess.h"
 #include "ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPProcess.h"
+#include "ProcessLib/TwoPhaseFlowWithPrho/CreateTwoPhaseFlowWithPrhoProcess.h"
 
 namespace detail
 {
@@ -442,6 +443,15 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config,
                     process_config, _curves);
         }
 
+        else if (type == "TWOPHASE_FLOW_PRHO")
+        {
+            process = ProcessLib::TwoPhaseFlowWithPrho::
+                createTwoPhaseFlowWithPrhoProcess(
+                    *_mesh_vec[0], std::move(jacobian_assembler),
+                    _process_variables, _parameters, integration_order,
+                    process_config, _curves);
+        }
+
         else
         {
             OGS_FATAL("Unknown process type: %s", type.c_str());
@@ -526,8 +536,8 @@ void ProjectData::parseCurves(
         BaseLib::insertIfKeyUniqueElseError(
             _curves,
             name,
-            MathLib::createPiecewiseLinearCurve
-                                 <MathLib::PiecewiseLinearInterpolation>(conf),
+            MathLib::createPiecewiseLinearCurve<
+                MathLib::PiecewiseLinearInterpolation>(conf),
             "The curve name is not unique.");
     }
 }

Documentation/ProjectFile/material/porous_medium/capillary_pressure/BrookCorey/t_sg_r.md

@@ -0,0 +1 @@
+A tag for the residual saturation in the gas phase.

Documentation/ProjectFile/material/porous_medium/capillary_pressure/vanGenuchten/t_has_regularized.md

@@ -0,0 +1 @@
+A tag for using regularized van Genuchten model.

Documentation/ProjectFile/material/porous_medium/capillary_pressure/vanGenuchten/t_sg_r.md

@@ -0,0 +1 @@
+A tag for the residual saturation of the nonwetting phase.

Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/c_TWOPHASE_FLOW_PRHO.md

@@ -0,0 +1 @@
+A tag for compositional based two-phase flow process, using liquid pressure and total mass density as process variables.

Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/i_material_property.md

@@ -0,0 +1 @@
+A tag for liquid and gas material properties of the two-phase flow process.

Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/porous_medium/i_porous_medium.md

@@ -0,0 +1 @@
+A tag for the properties of porous media.

Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/porous_medium/porous_medium/a_id.md

@@ -0,0 +1 @@
+It specifies the material ID of a porous medium. The material IDs of elements should be given in the mesh's cell data array "MaterialIDs".

Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/porous_medium/porous_medium/i_porous_medium.md

@@ -0,0 +1 @@
+A tag for the properties of a porous medium with material ID.

Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/porous_medium/porous_medium/relative_permeability/i_relative_permeability.md

@@ -0,0 +1 @@
+A tag for the relative permeability model, which holds two ids: one is for the Nonwetting phase and the other is for the Wetting phase.

Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/porous_medium/porous_medium/relative_permeability/relative_permeability/a_id.md

@@ -0,0 +1 @@
+A tag for the relative permeability model id: 0 is Nonwetting phase and 1 is Wetting phase.

Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/porous_medium/porous_medium/relative_permeability/relative_permeability/i_relative_permeability.md

@@ -0,0 +1 @@
+Defines the relative permeability model.

Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/porous_medium/porous_medium/t_capillary_pressure.md

@@ -0,0 +1 @@
+Capillary pressure model.

Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/porous_medium/porous_medium/t_permeability.md

@@ -0,0 +1 @@
+A tag for the intrinsic permeability model.

Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/porous_medium/porous_medium/t_porosity.md

@@ -0,0 +1 @@
+A tag for the porosity model.

Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/porous_medium/porous_medium/t_storage.md

@@ -0,0 +1 @@
+A tag for the storage model.

Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/t_fluid.md

@@ -0,0 +1 @@
+Defines fluid properties of gas and liquid phases.

Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/t_gas_density.md

@@ -0,0 +1 @@
+Density model of the gas phase.

Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/t_gas_viscosity.md

@@ -0,0 +1 @@
+Viscosity model of the gas phase.

Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/t_liquid_density.md

@@ -0,0 +1 @@
+Density model of the liquid phase.

Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/t_liquid_viscosity.md

@@ -0,0 +1 @@
+Viscosity model of the liquid phase.

Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/process_variables/i_process_variables.md

@@ -0,0 +1 @@
+The process variables for liquid pressure and overall mass density.

Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/process_variables/t_liquid_pressure.md

@@ -0,0 +1 @@
+Pressure for the liquid phase.

Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/process_variables/t_overall_mass_density.md

@@ -0,0 +1 @@
+Overall mass density of the light component

Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/t_diffusion_coeff_component_a.md

@@ -0,0 +1 @@
+Diffusion coefficient for the heavy component.

Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/t_diffusion_coeff_component_b.md

@@ -0,0 +1 @@
+Diffusion coefficient for the light component.

Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/t_mass_lumping.md

@@ -0,0 +1 @@
+A tag for using mass-lumping technique.

Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/t_specific_body_force.md

@@ -0,0 +1 @@
+\copydoc ProcessLib::TwoPhaseFlowWithPrho::TwoPhaseFlowWithPrhoProcessData::_specific_body_force

Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/t_temperature.md

@@ -0,0 +1 @@
+Constant temperature (K) for the simulation.

MaterialLib/PhysicalConstant.h

@@ -74,6 +74,13 @@ const double O2 = 0.032;  ///< kg \per{mol}
  * - Source: http://www.engineeringtoolbox.com/molecular-mass-air-d_679.html
  */
 const double Air = 0.02897;  ///< kg \per{mol}
+
+/**
+ * Hydrogen.
+ *
+ * - Source: https://pubchem.ncbi.nlm.nih.gov/compound/Hydrogen
+ */
+const double H2 = 0.002016;  ///< kg \per{mol}
 }  // namespace MolarMass
 
 /**
@@ -86,5 +93,18 @@ namespace SpecificGasConstant
 /// Specific gas constant for water vapour.
 const double WaterVapour = IdealGasConstant / MolarMass::Water;  // = 461.504;
 }  // namespace SpecificGasConstant
+/**
+* Henry's law constant
+* It is defined here as the ratio of the aqueous-phase concentration of a
+* chemical to its equilibrium partial pressure in the gas phase.
+*/
+namespace HenryConstant
+{
+/**Henry constant for hydrogen, Please refer to
+* De Nevers N. Physical and chemical equilibrium for chemical engineers[M].
+* John Wiley & Sons, 2012.
+*/
+double const HenryConstantH2 = 7.65e-6;  /// mol/Pa./m3
+}  // namespace HenryConstant
 }  // namespace PhysicalConstant
 }  // namespace MaterialLib

MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/BrookCoreyCapillaryPressureSaturation.h

@@ -44,14 +44,15 @@ class BrookCoreyCapillaryPressureSaturation final
     /**
      * @param pb     Entry pressure, \f$ p_b \f$
      * @param Sr     Residual saturation, \f$ S_r \f$
+     * @param Sg_r     Residual saturation of gas phase, \f$ Sg_r \f$
      * @param Smax   Maximum saturation, \f$ S_{\mbox{max}} \f$
      * @param m      Exponent, \f$ m \f$
      * @param Pc_max Maximum capillary pressure, \f$ P_c^{\mbox{max}}\f$
      */
     BrookCoreyCapillaryPressureSaturation(const double pb, const double Sr,
-                                          const double Smax, const double m,
-                                          const double Pc_max)
-        : CapillaryPressureSaturation(Sr, Smax, Pc_max), _pb(pb), _m(m)
+                                          const double Sg_r, const double Smax,
+                                          const double m, const double Pc_max)
+        : CapillaryPressureSaturation(Sr, Sg_r, Smax, Pc_max), _pb(pb), _m(m)
     {
     }
 

MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CapillaryPressureSaturation.h

@@ -23,12 +23,16 @@ class CapillaryPressureSaturation
 public:
     /**
      * @param Sr     Residual saturation, \f$ S_r \f$
+     * @param Sg_r     Residual saturation of gas phase, \f$ Sg_r \f$
      * @param Smax   Maximum saturation, \f$ S_{\mbox{max}} \f$
      * @param Pc_max Maximum capillary pressure, \f$ P_c^{\mbox{max}}\f$
      */
-    CapillaryPressureSaturation(const double Sr, const double Smax,
-                                const double Pc_max)
-        : _saturation_r(Sr), _saturation_max(Smax), _pc_max(Pc_max)
+    CapillaryPressureSaturation(const double Sr, const double Sg_r,
+                                const double Smax, const double Pc_max)
+        : _saturation_r(Sr),
+          _saturation_nonwet_r(Sg_r),
+          _saturation_max(Smax),
+          _pc_max(Pc_max)
     {
     }
     virtual ~CapillaryPressureSaturation() = default;
@@ -38,17 +42,18 @@ class CapillaryPressureSaturation
 
     /// Get capillary pressure.
     virtual double getCapillaryPressure(const double saturation) const = 0;
-
     /// Get saturation.
     virtual double getSaturation(const double capillary_ressure) const = 0;
 
     /// Get the derivative of the capillary pressure with respect to saturation
     virtual double getdPcdS(const double saturation) const = 0;
 
 protected:
-    const double _saturation_r;    ///< Residual saturation.
-    const double _saturation_max;  ///< Maximum saturation.
-    const double _pc_max;          ///< Maximum capillaray pressure
+    const double _saturation_r;         ///< Residual saturation.
+    const double _saturation_nonwet_r;  ///< Residual saturation of nonwetting
+                                        ///phase (optional).
+    const double _saturation_max;       ///< Maximum saturation.
+    const double _pc_max;               ///< Maximum capillaray pressure
 
     /** A small number for an offset:
      *  1. to set the bound of S, the saturation, such that

MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CapillaryPressureSaturationCurve.h

@@ -31,6 +31,7 @@ class CapillaryPressureSaturationCurve final
         std::unique_ptr<MathLib::PiecewiseLinearMonotonicCurve>&& curve_data)
         : CapillaryPressureSaturation(
               curve_data->getSupportMin(),
+              1.0 - curve_data->getSupportMax(),
               curve_data->getSupportMax(),
               curve_data->getValue(curve_data->getSupportMin())),
           _curve_data(std::move(curve_data))

MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.cpp

@@ -44,6 +44,16 @@ static std::unique_ptr<CapillaryPressureSaturation> createBrookCorey(
     //! \ogs_file_param{material__porous_medium__capillary_pressure__BrookCorey__sr}
     const double Sr = config.getConfigParameter<double>("sr");
 
+    double Sg_r = 0.0;
+    //! \ogs_file_param{material__porous_medium__capillary_pressure__BrookCorey__sg_r}
+    if (auto const Sg_r_ptr = config.getConfigParameterOptional<double>("sg_r"))
+    {
+        DBUG(
+            "Using value %g for nonwetting phase residual saturation in "
+            "capillary pressure model.",
+            (*Sg_r_ptr));
+        Sg_r = *Sg_r_ptr;
+    }
     //! \ogs_file_param{material__porous_medium__capillary_pressure__BrookCorey__smax}
     const double Smax = config.getConfigParameter<double>("smax");
 
@@ -59,7 +69,8 @@ static std::unique_ptr<CapillaryPressureSaturation> createBrookCorey(
     const double Pc_max = config.getConfigParameter<double>("pc_max");
 
     return std::unique_ptr<CapillaryPressureSaturation>(
-        new BrookCoreyCapillaryPressureSaturation(pd, Sr, Smax, m, Pc_max));
+        new BrookCoreyCapillaryPressureSaturation(
+            pd, Sr, Sg_r, Smax, m, Pc_max));
 }
 
 /**
@@ -79,6 +90,17 @@ static std::unique_ptr<CapillaryPressureSaturation> createVanGenuchten(
     //! \ogs_file_param{material__porous_medium__capillary_pressure__vanGenuchten__sr}
     const double Sr = config.getConfigParameter<double>("sr");
 
+    double Sg_r = 0.0;
+    //! \ogs_file_param{material__porous_medium__capillary_pressure__vanGenuchten__sg_r}
+    if (auto const Sg_r_ptr = config.getConfigParameterOptional<double>("sg_r"))
+    {
+        DBUG(
+            "Using value %g for nonwetting phase residual saturation in "
+            "capillary pressure model.",
+            (*Sg_r_ptr));
+        Sg_r = *Sg_r_ptr;
+    }
+
     //! \ogs_file_param{material__porous_medium__capillary_pressure__vanGenuchten__smax}
     const double Smax = config.getConfigParameter<double>("smax");
 
@@ -93,8 +115,18 @@ static std::unique_ptr<CapillaryPressureSaturation> createVanGenuchten(
     //! \ogs_file_param{material__porous_medium__capillary_pressure__vanGenuchten__pc_max}
     const double Pc_max = config.getConfigParameter<double>("pc_max");
 
+    bool has_regularized = false;
+    if (auto const has_regularized_conf =
+            //! \ogs_file_param{material__porous_medium__capillary_pressure__vanGenuchten__has_regularized}
+        config.getConfigParameterOptional<bool>("has_regularized"))
+    {
+        DBUG("capillary pressure model: %s",
+             (*has_regularized_conf) ? "true" : "false");
+        has_regularized = *has_regularized_conf;
+    }
     return std::unique_ptr<CapillaryPressureSaturation>(
-        new VanGenuchtenCapillaryPressureSaturation(pd, Sr, Smax, m, Pc_max));
+        new VanGenuchtenCapillaryPressureSaturation(
+            pd, Sr, Sg_r, Smax, m, Pc_max, has_regularized));
 }
 
 std::unique_ptr<CapillaryPressureSaturation> createCapillaryPressureModel(

MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/VanGenuchtenCapillaryPressureSaturation.cpp

@@ -23,6 +23,25 @@ namespace PorousMedium
 double VanGenuchtenCapillaryPressureSaturation::getCapillaryPressure(
     const double saturation) const
 {
+    if (_has_regularized)
+    {
+        double Sg = 1 - saturation;
+        if (Sg <= 1 - _saturation_r && Sg >= _saturation_nonwet_r)
+        {
+            return getPcBarvGSg(Sg);
+        }
+        else if (Sg < _saturation_nonwet_r)
+        {
+            return getPcBarvGSg(_saturation_nonwet_r) +
+                   getdPcdSvGBar(_saturation_nonwet_r) *
+                       (Sg - _saturation_nonwet_r);
+        }
+        else
+        {
+            return getPcBarvGSg(1 - _saturation_r) +
+                   getdPcdSvGBar(1 - _saturation_r) * (Sg - 1 + _saturation_r);
+        }
+    }
     const double S =
         MathLib::limitValueInInterval(saturation, _saturation_r + _minor_offset,
                                       _saturation_max - _minor_offset);
@@ -46,6 +65,22 @@ double VanGenuchtenCapillaryPressureSaturation::getSaturation(
 double VanGenuchtenCapillaryPressureSaturation::getdPcdS(
     const double saturation) const
 {
+    if (_has_regularized)
+    {
+        double const Sg = 1 - saturation;
+        if (Sg >= _saturation_nonwet_r && Sg <= 1 - _saturation_r)
+        {
+            return -getdPcdSvGBar(Sg);
+        }
+        else if (Sg < _saturation_nonwet_r)
+        {
+            return -getdPcdSvGBar(_saturation_nonwet_r);
+        }
+        else
+        {
+            return -getdPcdSvGBar(1 - _saturation_r);
+        }
+    }
     const double S =
         MathLib::limitValueInInterval(saturation, _saturation_r + _minor_offset,
                                       _saturation_max - _minor_offset);
@@ -54,6 +89,50 @@ double VanGenuchtenCapillaryPressureSaturation::getdPcdS(
     const double val2 = std::pow(val1 - 1.0, -_m);
     return _pb * (_m - 1.0) * val1 * val2 / (_m * (S - _saturation_r));
 }
+/// Regularized van Genuchten capillary pressure-saturation Model
+double VanGenuchtenCapillaryPressureSaturation::getPcBarvGSg(double Sg) const
+{
+    double const Sg_r = CapillaryPressureSaturation::_saturation_nonwet_r;
+    double const S_lr = CapillaryPressureSaturation::_saturation_r;
+    double const S_bar = getSBar(Sg);
+    return getPcvGSg(S_bar) - getPcvGSg(Sg_r + (1 - Sg_r - S_lr) * _xi / 2);
+}
+/// Regularized van Genuchten capillary pressure-saturation Model
+double VanGenuchtenCapillaryPressureSaturation::getSBar(double Sg) const
+{
+    double const Sg_r = CapillaryPressureSaturation::_saturation_nonwet_r;
+    double const S_lr = CapillaryPressureSaturation::_saturation_r;
+    return Sg_r + (1 - _xi) * (Sg - Sg_r) + 0.5 * _xi * (1 - Sg_r - S_lr);
+}
+///  van Genuchten capillary pressure-saturation Model
+double VanGenuchtenCapillaryPressureSaturation::getPcvGSg(double Sg) const
+{
+    double const Sg_r = CapillaryPressureSaturation::_saturation_nonwet_r;
+    double const S_lr = CapillaryPressureSaturation::_saturation_r;
+    double const S_le = (1 - Sg - S_lr) /
+                        (1 - Sg_r - CapillaryPressureSaturation::_saturation_r);
+    return _pb * std::pow(std::pow(S_le, (-1.0 / _m)) - 1.0, 1.0 - _m);
+}
+/// derivative dPCdS based on regularized van Genuchten capillary
+/// pressure-saturation Model
+double VanGenuchtenCapillaryPressureSaturation::getdPcdSvGBar(double Sg) const
+{
+    double S_bar = getSBar(Sg);
+    return getdPcdSvG(S_bar) * (1 - _xi);
+}
+/// derivative dPCdS based on standard van Genuchten capillary
+/// pressure-saturation Model
+double VanGenuchtenCapillaryPressureSaturation::getdPcdSvG(
+    const double Sg) const
+{
+    double const Sg_r = CapillaryPressureSaturation::_saturation_nonwet_r;
+    double const S_lr = CapillaryPressureSaturation::_saturation_r;
+    double const nn = 1 / (1 - _m);
+    double const S_le = (1 - Sg - S_lr) / (1 - Sg_r - S_lr);
+    return _pb * (1 / (_m * nn)) * (1 / (1 - S_lr - Sg_r)) *
+           std::pow(std::pow(S_le, (-1 / _m)) - 1, (1 / nn) - 1) *
+           std::pow(S_le, (-1 / _m)) / S_le;
+}
 
 }  // end namespace
 }  // end namespace