Revision deea81a23b76b9fc96a0e411a8e72670ca85eefc authored by Dmitry Yu. Naumov on 29 July 2021, 10:55:10 UTC, committed by Dmitry Yu. Naumov on 29 July 2021, 10:55:10 UTC
Enabled integration point data I/O in THM See merge request ogs/ogs!3718
TwoPhaseFlowWithPPMaterialProperties.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
*/
#include "TwoPhaseFlowWithPPMaterialProperties.h"
#include <boost/math/special_functions/pow.hpp>
#include <utility>
#include "BaseLib/Logging.h"
#include "MaterialLib/Fluid/FluidProperty.h"
#include "MaterialLib/PorousMedium/Porosity/Porosity.h"
#include "MaterialLib/PorousMedium/Storage/Storage.h"
#include "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CapillaryPressureSaturation.h"
#include "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.h"
#include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h"
#include "MeshLib/Mesh.h"
#include "MeshLib/PropertyVector.h"
#include "ParameterLib/Parameter.h"
#include "ParameterLib/SpatialPosition.h"
namespace MaterialLib
{
namespace TwoPhaseFlowWithPP
{
TwoPhaseFlowWithPPMaterialProperties::TwoPhaseFlowWithPPMaterialProperties(
MeshLib::PropertyVector<int> const* material_ids,
std::unique_ptr<MaterialLib::Fluid::FluidProperty>&& liquid_density,
std::unique_ptr<MaterialLib::Fluid::FluidProperty>&& liquid_viscosity,
std::unique_ptr<MaterialLib::Fluid::FluidProperty>&& gas_density,
std::unique_ptr<MaterialLib::Fluid::FluidProperty>&& gas_viscosity,
std::vector<std::unique_ptr<MaterialLib::PorousMedium::Permeability>>&&
intrinsic_permeability_models,
std::vector<std::unique_ptr<MaterialLib::PorousMedium::Porosity>>&&
porosity_models,
std::vector<std::unique_ptr<MaterialLib::PorousMedium::Storage>>&&
storage_models,
std::vector<std::unique_ptr<
MaterialLib::PorousMedium::CapillaryPressureSaturation>>&&
capillary_pressure_models,
std::vector<
std::unique_ptr<MaterialLib::PorousMedium::RelativePermeability>>&&
relative_permeability_models)
: _material_ids(material_ids),
_liquid_density(std::move(liquid_density)),
_liquid_viscosity(std::move(liquid_viscosity)),
_gas_density(std::move(gas_density)),
_gas_viscosity(std::move(gas_viscosity)),
_intrinsic_permeability_models(std::move(intrinsic_permeability_models)),
_porosity_models(std::move(porosity_models)),
_storage_models(std::move(storage_models)),
_capillary_pressure_models(std::move(capillary_pressure_models)),
_relative_permeability_models(std::move(relative_permeability_models))
{
DBUG("Create material properties for Two-Phase flow with PP model.");
}
int TwoPhaseFlowWithPPMaterialProperties::getMaterialID(
const std::size_t element_id) const
{
if (_material_ids == nullptr || _material_ids->empty())
{
return 0;
}
assert(element_id < _material_ids->size());
return (*_material_ids)[element_id];
}
double TwoPhaseFlowWithPPMaterialProperties::getLiquidDensity(
const double p, const double T) const
{
ArrayType vars;
vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)] = T;
vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::p)] = p;
return _liquid_density->getValue(vars);
}
double TwoPhaseFlowWithPPMaterialProperties::getGasDensity(const double p,
const double T) const
{
ArrayType vars;
vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)] = T;
vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::p)] = p;
return _gas_density->getValue(vars);
}
double TwoPhaseFlowWithPPMaterialProperties::getLiquidViscosity(
const double p, const double T) const
{
ArrayType vars;
vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)] = T;
vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::p)] = p;
return _liquid_viscosity->getValue(vars);
}
double TwoPhaseFlowWithPPMaterialProperties::getGasViscosity(
const double p, const double T) const
{
ArrayType vars;
vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)] = T;
vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::p)] = p;
return _gas_viscosity->getValue(vars);
}
Eigen::MatrixXd TwoPhaseFlowWithPPMaterialProperties::getPermeability(
const int material_id, const double t,
const ParameterLib::SpatialPosition& pos, const int /*dim*/) const
{
return _intrinsic_permeability_models[material_id]->getValue(t, pos, 0, 0);
}
double TwoPhaseFlowWithPPMaterialProperties::getPorosity(
const int material_id, const double t,
const ParameterLib::SpatialPosition& pos, const double /*p*/,
const double T, const double porosity_variable) const
{
return _porosity_models[material_id]->getValue(t, pos, porosity_variable,
T);
}
double TwoPhaseFlowWithPPMaterialProperties::getSaturation(
const int material_id, const double /*t*/,
const ParameterLib::SpatialPosition& /*pos*/, const double /*p*/,
const double /*T*/, const double pc) const
{
return _capillary_pressure_models[material_id]->getSaturation(pc);
}
double TwoPhaseFlowWithPPMaterialProperties::getCapillaryPressure(
const int material_id, const double /*t*/,
const ParameterLib::SpatialPosition& /*pos*/, const double /*p*/,
const double /*T*/, const double saturation) const
{
return _capillary_pressure_models[material_id]->getCapillaryPressure(
saturation);
}
double TwoPhaseFlowWithPPMaterialProperties::getSaturationDerivative(
const int material_id, const double /*t*/,
const ParameterLib::SpatialPosition& /*pos*/, const double /*p*/,
const double /*T*/, const double saturation) const
{
const double dpcdsw =
_capillary_pressure_models[material_id]->getdPcdS(saturation);
return 1 / dpcdsw;
}
double TwoPhaseFlowWithPPMaterialProperties::getNonwetRelativePermeability(
const double /*t*/, const ParameterLib::SpatialPosition& /*pos*/,
const double /*p*/, const double /*T*/, const double saturation) const
{
if (saturation < 0.)
{
return 1.0;
}
if (saturation > 1)
{
return 0.0;
}
return boost::math::pow<3>(1 - saturation);
}
double TwoPhaseFlowWithPPMaterialProperties::getWetRelativePermeability(
const double /*t*/, const ParameterLib::SpatialPosition& /*pos*/,
const double /*p*/, const double /*T*/, const double saturation) const
{
if (saturation < 0)
{
return 0.0;
}
if (saturation > 1)
{
return 1.0;
}
return boost::math::pow<3>(saturation);
}
} // namespace TwoPhaseFlowWithPP
} // namespace MaterialLib
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...