Commit 385d7007 authored by Dmitry Yu. Naumov's avatar Dmitry Yu. Naumov
Browse files

Merge branch 'MPL_Permmodel_Tutorial_TEC' into 'master'

Gas Pressure Permeability Model

See merge request ogs/ogs!3439
parents 74e9377e ac4f8c5c
\copydoc MaterialPropertyLib::GasPressureDependentPermeability
\ No newline at end of file
\copydoc MaterialPropertyLib::GasPressureDependentPermeability::a1_
\copydoc MaterialPropertyLib::GasPressureDependentPermeability::a2_
\copydoc MaterialPropertyLib::GasPressureDependentPermeability::pressure_threshold_
......@@ -87,6 +87,12 @@ std::unique_ptr<MaterialPropertyLib::Property> createProperty(
geometry_dimension, config, parameters, local_coordinate_system);
}
if (boost::iequals(property_type, "GasPressureDependentPermeability"))
{
return createGasPressureDependentPermeability(
geometry_dimension, config, parameters, local_coordinate_system);
}
if (boost::iequals(property_type, "EmbeddedFracturePermeability"))
{
return createEmbeddedFracturePermeability(geometry_dimension, config);
......
/**
* \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 January 12, 2021, x:xx AM
*/
#include "CreateGasPressureDependentPermeability.h"
#include <string>
#include "BaseLib/ConfigTree.h"
#include "BaseLib/Error.h"
#include "GasPressureDependentPermeability.h"
#include "Parameter.h"
#include "ParameterLib/CoordinateSystem.h"
#include "ParameterLib/Parameter.h"
#include "ParameterLib/Utils.h"
namespace MaterialPropertyLib
{
std::unique_ptr<Property> createGasPressureDependentPermeability(
int const geometry_dimension,
BaseLib::ConfigTree const& config,
std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
ParameterLib::CoordinateSystem const* const local_coordinate_system)
{
if ((geometry_dimension != 2) && (geometry_dimension != 3))
{
OGS_FATAL(
"The GasPressureDependentPermeability is implemented only "
"for 2D or 3D problems");
}
//! \ogs_file_param{properties__property__type}
config.checkConfigParameter("type", "GasPressureDependentPermeability");
// Second access for storage.
//! \ogs_file_param{properties__property__name}
auto property_name = config.peekConfigParameter<std::string>("name");
DBUG("Create GasPressureDependentPermeability property {:s}.",
property_name);
std::string const& parameter_name =
//! \ogs_file_param{properties__property__GasPressureDependentPermeability__initial_permeability}
config.getConfigParameter<std::string>("initial_permeability");
auto const& parameter_k0 = ParameterLib::findParameter<double>(
parameter_name, parameters, 0, nullptr);
auto const a1 =
//! \ogs_file_param{properties__property__GasPressureDependentPermeability__a1}
config.getConfigParameter<double>("a1");
auto const a2 =
//! \ogs_file_param{properties__property__GasPressureDependentPermeability__a2}
config.getConfigParameter<double>("a2");
auto const pressure_threshold =
//! \ogs_file_param{properties__property__GasPressureDependentPermeability__pressure_threshold}
config.getConfigParameter<double>("pressure_threshold");
auto const minimum_permeability =
//! \ogs_file_param{properties__property__GasPressureDependentPermeability__minimum_permeability}
config.getConfigParameter<double>("minimum_permeability");
auto const maximum_permeability =
//! \ogs_file_param{properties__property__GasPressureDependentPermeability__maximum_permeability}
config.getConfigParameter<double>("maximum_permeability");
if (minimum_permeability > maximum_permeability)
{
OGS_FATAL(
"The value of minimum_permeability of {:e} is larger that the "
"value of maximum_permeability of {:e} in "
"GasPressureDependentPermeability",
minimum_permeability, maximum_permeability);
}
if (geometry_dimension == 2)
{
return std::make_unique<GasPressureDependentPermeability<2>>(
std::move(property_name), parameter_k0, a1, a2, pressure_threshold,
minimum_permeability, maximum_permeability,
local_coordinate_system);
}
return std::make_unique<GasPressureDependentPermeability<3>>(
std::move(property_name), parameter_k0, a1, a2, pressure_threshold,
minimum_permeability, maximum_permeability, local_coordinate_system);
}
} // namespace MaterialPropertyLib
/**
* \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 January 12, 2021, x:xx AM
*/
#pragma once
#include <memory>
#include <vector>
namespace BaseLib
{
class ConfigTree;
}
namespace ParameterLib
{
struct CoordinateSystem;
struct ParameterBase;
} // namespace ParameterLib
namespace MaterialPropertyLib
{
class Property;
}
namespace MaterialPropertyLib
{
std::unique_ptr<Property> createGasPressureDependentPermeability(
int const geometry_dimension,
BaseLib::ConfigTree const& config,
std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
ParameterLib::CoordinateSystem const* const local_coordinate_system);
} // namespace MaterialPropertyLib
......@@ -25,6 +25,7 @@
#include "CreateDupuitPermeability.h"
#include "CreateEmbeddedFracturePermeability.h"
#include "CreateExponential.h"
#include "CreateGasPressureDependentPermeability.h"
#include "CreateIdealGasLaw.h"
#include "CreateKozenyCarmanModel.h"
#include "CreateLinear.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 January 12, 2021, x:xx AM
*/
#include "GasPressureDependentPermeability.h"
#include <cmath>
#include <limits>
#include "BaseLib/Error.h"
#include "MaterialLib/MPL/Medium.h"
#include "MathLib/MathTools.h"
#include "ParameterLib/CoordinateSystem.h"
#include "ParameterLib/Parameter.h"
namespace MaterialPropertyLib
{
template <int DisplacementDim>
GasPressureDependentPermeability<DisplacementDim>::
GasPressureDependentPermeability(
std::string name, ParameterLib::Parameter<double> const& k0,
double const a1, double const a2, double const pressure_threshold,
double const minimum_permeability, double const maximum_permeability,
ParameterLib::CoordinateSystem const* const local_coordinate_system)
: k0_(k0),
a1_(a1),
a2_(a2),
pressure_threshold_(pressure_threshold),
minimum_permeability_(minimum_permeability),
maximum_permeability_(maximum_permeability),
local_coordinate_system_(local_coordinate_system)
{
name_ = std::move(name);
}
template <int DisplacementDim>
void GasPressureDependentPermeability<DisplacementDim>::checkScale() const
{
if (!std::holds_alternative<Medium*>(scale_))
{
OGS_FATAL(
"The property 'GasPressureDependentPermeability' is "
"implemented on the 'medium' scale only.");
}
}
template <int DisplacementDim>
PropertyDataType GasPressureDependentPermeability<DisplacementDim>::value(
VariableArray const& variable_array,
ParameterLib::SpatialPosition const& pos, double const t,
double const /*dt*/) const
{
double const gas_pressure = std::get<double>(
variable_array[static_cast<int>(Variable::phase_pressure)]);
auto k_data = k0_(t, pos);
double const factor = (gas_pressure <= pressure_threshold_)
? (1.0 + a1_ * gas_pressure)
: (a2_ * (gas_pressure - pressure_threshold_) +
1.0 + a1_ * pressure_threshold_);
for (auto& k_i : k_data)
{
k_i = std::clamp(k_i * factor, minimum_permeability_,
maximum_permeability_);
}
// Local coordinate transformation is only applied for the case that the
// initial intrinsic permeability is given with orthotropic assumption.
if (local_coordinate_system_ && (k_data.size() == DisplacementDim))
{
Eigen::Matrix<double, DisplacementDim, DisplacementDim> const e =
local_coordinate_system_->transformation<DisplacementDim>(pos);
Eigen::Matrix<double, DisplacementDim, DisplacementDim> k =
Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
for (int i = 0; i < DisplacementDim; ++i)
{
Eigen::Matrix<double, DisplacementDim, DisplacementDim> const
ei_otimes_ei = e.col(i) * e.col(i).transpose();
k += k_data[i] * ei_otimes_ei;
}
return k;
}
return fromVector(k_data);
}
template <int DisplacementDim>
PropertyDataType GasPressureDependentPermeability<DisplacementDim>::dValue(
VariableArray const& /*variable_array*/, Variable const /*variable*/,
ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
double const /*dt*/) const
{
OGS_FATAL(
"The derivative of the intrinsic permeability of "
"GasPressureDependentPermeability"
"is not implemented.");
}
template class GasPressureDependentPermeability<2>;
template class GasPressureDependentPermeability<3>;
} // namespace MaterialPropertyLib
/**
* \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 January 12, 2021, x:xx AM
*/
#pragma once
#include "MaterialLib/MPL/Property.h"
#include "MaterialLib/MPL/VariableType.h"
namespace ParameterLib
{
struct CoordinateSystem;
template <typename T>
struct Parameter;
} // namespace ParameterLib
namespace MaterialPropertyLib
{
/**
* \brief A gas pressure dependent intrinsic permeability model.
*
* The model was proposed
* in \cite xu2011simulation and it was further investigated
* in \cite xu2013coupled .
*
* The model takes the form of
* \f[\mathbf{k} = f(p_g)\, \mathbf{k}_0
* =
* \left\{
* \begin{array}{cc}
* (1+a_1 p_g) \mathbf{k}_0 & p_g \leq p_\text{thr} \\
* (a_2 (p_g - p_\text{thr}) + 1 + a_1 p_\text{thr}) \mathbf{k}_0 & p_g >
* p_\text{thr} \end{array} \right.\f]
* where \f$\mathbf{k}\f$ is the permeability, \f$\mathbf{k}_0\f$ is the
* initial intrinsic permeability, \f$p_g\f$ is the gas pressure, \f$a_1\f$,
* \f$a_2\f$ and \f$p_\text{thr}\f$ are three parameters.
*
*/
template <int DisplacementDim>
class GasPressureDependentPermeability final : public Property
{
public:
GasPressureDependentPermeability(
std::string name, ParameterLib::Parameter<double> const& k0,
double const a1, double const a2, double const pressure_threshold,
double const minimum_permeability, double const maximum_permeability,
ParameterLib::CoordinateSystem const* const local_coordinate_system);
void checkScale() const override;
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:
/// Initial intrinsic permeability.
ParameterLib::Parameter<double> const& k0_;
/// A factor taking into account the influence of the gas pressure on the
/// permeability.
double const a1_;
/// A factor taking into account the influence of the gas overpressure
/// (difference of the gas pressure and the threshold pressure) on the
/// permeability.
double const a2_;
/// The threshold pressure which defines the low and high gas pressure
/// domain. Different permeability functions are defined for each domain.
double const pressure_threshold_;
double const minimum_permeability_;
double const maximum_permeability_;
ParameterLib::CoordinateSystem const* const local_coordinate_system_;
};
extern template class GasPressureDependentPermeability<2>;
extern template class GasPressureDependentPermeability<3>;
} // namespace MaterialPropertyLib
......@@ -23,6 +23,7 @@
#include "DupuitPermeability.h"
#include "EmbeddedFracturePermeability.h"
#include "Exponential.h"
#include "GasPressureDependentPermeability.h"
#include "IdealGasLaw.h"
#include "Linear.h"
#include "Parameter.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 January 12, 2021, x:xx AM
*/
#include <gtest/gtest.h>
#include <Eigen/Eigen>
#include <boost/math/constants/constants.hpp>
#include "MaterialLib/MPL/Medium.h"
#include "MaterialLib/MPL/Properties/CreateGasPressureDependentPermeability.h"
#include "MaterialLib/MPL/Properties/GasPressureDependentPermeability.h"
#include "MaterialLib/MPL/Utils/FormEigenTensor.h"
#include "ParameterLib/ConstantParameter.h"
#include "TestMPL.h"
#include "Tests/TestTools.h"
TEST(MaterialPropertyLib, GasPressureDependentPermeability)
{
ParameterLib::ConstantParameter<double> const k0("k0", 1.e-20);
double const a1 = 0.125;
double const a2 = 152.0;
double const pressure_threshold = 3.2;
double const min_permeability = 1.e-22;
double const max_permeability = 1.e-10;
auto const k_model = MPL::GasPressureDependentPermeability<3>(
"k_gas", k0, a1, a2, pressure_threshold, min_permeability,
max_permeability, nullptr);
ParameterLib::SpatialPosition const pos;
double const t = std::numeric_limits<double>::quiet_NaN();
double const dt = std::numeric_limits<double>::quiet_NaN();
MPL::VariableArray vars;
/// For gas pressure smaller than threshold value.
{
double const p_gas = 2.5;
vars[static_cast<int>(MPL::Variable::phase_pressure)] = p_gas;
auto const k = MPL::formEigenTensor<3>(k_model.value(vars, pos, t, dt));
double const k_expected = 1.312500000000000000e-20;
ASSERT_LE(std::fabs(k_expected - k(0, 0)) / k_expected, 1e-10)
<< "for expected permeability with gas pressure below threshold"
<< k_expected
<< " and for computed permeability with gas pressure below "
"threshold "
<< k(0, 0);
}
/// For gas pressure bigger than threshold value.
{
double const p_gas = 4.5;
vars[static_cast<int>(MPL::Variable::phase_pressure)] = p_gas;
auto const k = MPL::formEigenTensor<3>(k_model.value(vars, pos, t, dt));
double const k_expected = 1.990000000000000000000e-18;
ASSERT_LE(std::fabs(k_expected - k(0, 0)) / k_expected, 1e-10)
<< "for expected permeability with gas pressure above threshold "
<< k_expected
<< " and for computed permeability with gas pressure above "
"threshold "
<< k(0, 0);
}
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment