From 4225d1c5a61f9a271abdc50c245d617c96aa2afd Mon Sep 17 00:00:00 2001
From: Florian Zill <florian.zill@ufz.de>
Date: Fri, 23 Aug 2019 15:28:10 +0200
Subject: [PATCH] [HM] specific storage is now calculated added specific gas
constant added Fluid compressibility added Fluid type added bbulk modulus
removed specific storage
---
.../t_fluid_compressibility.md | 1 +
.../process/HYDRO_MECHANICS/t_fluid_type.md | 1 +
.../t_specific_gas_constant.md | 1 +
.../CreateHydroMechanicsProcess.cpp | 49 +++++++++++-----
.../HydroMechanics/HydroMechanicsFEM-impl.h | 40 ++++++++++---
.../HydroMechanicsProcessData.h | 56 ++++++++++++++++++-
6 files changed, 122 insertions(+), 26 deletions(-)
create mode 100644 Documentation/ProjectFile/prj/processes/process/HYDRO_MECHANICS/t_fluid_compressibility.md
create mode 100644 Documentation/ProjectFile/prj/processes/process/HYDRO_MECHANICS/t_fluid_type.md
create mode 100644 Documentation/ProjectFile/prj/processes/process/HYDRO_MECHANICS/t_specific_gas_constant.md
diff --git a/Documentation/ProjectFile/prj/processes/process/HYDRO_MECHANICS/t_fluid_compressibility.md b/Documentation/ProjectFile/prj/processes/process/HYDRO_MECHANICS/t_fluid_compressibility.md
new file mode 100644
index 00000000000..8601db2de06
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/HYDRO_MECHANICS/t_fluid_compressibility.md
@@ -0,0 +1 @@
+\copydoc ProcessLib::HydroMechanics::HydroMechanicsProcessData::fluid_compressibility
diff --git a/Documentation/ProjectFile/prj/processes/process/HYDRO_MECHANICS/t_fluid_type.md b/Documentation/ProjectFile/prj/processes/process/HYDRO_MECHANICS/t_fluid_type.md
new file mode 100644
index 00000000000..24e00605c57
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/HYDRO_MECHANICS/t_fluid_type.md
@@ -0,0 +1 @@
+\copydoc ProcessLib::HydroMechanics::HydroMechanicsProcessData::fluid_type
diff --git a/Documentation/ProjectFile/prj/processes/process/HYDRO_MECHANICS/t_specific_gas_constant.md b/Documentation/ProjectFile/prj/processes/process/HYDRO_MECHANICS/t_specific_gas_constant.md
new file mode 100644
index 00000000000..45897fcdfa7
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/HYDRO_MECHANICS/t_specific_gas_constant.md
@@ -0,0 +1 @@
+\copydoc ProcessLib::HydroMechanics::HydroMechanicsProcessData::specific_gas_constant
diff --git a/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp b/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp
index 80b00d24701..e937977f35d 100644
--- a/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp
+++ b/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp
@@ -117,15 +117,6 @@ std::unique_ptr<Process> createHydroMechanicsProcess(
DBUG("Use '%s' as intrinsic conductivity parameter.",
intrinsic_permeability.name.c_str());
- // Storage coefficient
- auto& specific_storage = ParameterLib::findParameter<double>(
- config,
- //! \ogs_file_param_special{prj__processes__process__HYDRO_MECHANICS__specific_storage}
- "specific_storage", parameters, 1, &mesh);
-
- DBUG("Use '%s' as storage coefficient parameter.",
- specific_storage.name.c_str());
-
// Fluid viscosity
auto& fluid_viscosity = ParameterLib::findParameter<double>(
config,
@@ -183,18 +174,46 @@ std::unique_ptr<Process> createHydroMechanicsProcess(
}
// Reference temperature
- const auto& reference_temperature =
+ double const reference_temperature =
//! \ogs_file_param{prj__processes__process__HYDRO_MECHANICS__reference_temperature}
config.getConfigParameter<double>(
"reference_temperature", std::numeric_limits<double>::quiet_NaN());
+ DBUG("Use 'reference_temperature' as reference temperature.");
+
+ //Specific gas constant
+ double const specific_gas_constant =
+ //! \ogs_file_param{prj__processes__process__HYDRO_MECHANICS__specific_gas_constant}
+ config.getConfigParameter<double>(
+ "specific_gas_constant", std::numeric_limits<double>::quiet_NaN());
+ DBUG("Use 'specific_gas_constant' as specific gas constant.");
+
+ // Fluid compressibility
+ double const fluid_compressibility =
+ //! \ogs_file_param{prj__processes__process__HYDRO_MECHANICS__fluid_compressibility}
+ config.getConfigParameter<double>(
+ "fluid_compressibility", std::numeric_limits<double>::quiet_NaN());
+ DBUG("Use 'fluid_compressibility' as fluid compressibility parameter.");
+
+ auto const fluid_type =
+ //! \ogs_file_param{prj__processes__process__HYDRO_MECHANICS__fluid_type}
+ FluidType::strToFluidType(
+ config.getConfigParameter<std::string>("fluid_type"));
+ DBUG("Use 'fluid_type' as fluid type parameter.");
+
+ if (FluidType::checkRequiredParams(fluid_type, fluid_compressibility,
+ reference_temperature, specific_gas_constant) == false)
+ {
+ OGS_FATAL(FluidType::getErrorMsg(fluid_type));
+ }
HydroMechanicsProcessData<DisplacementDim> process_data{
materialIDs(mesh), std::move(solid_constitutive_relations),
- intrinsic_permeability, specific_storage,
- fluid_viscosity, fluid_density,
- biot_coefficient, porosity,
- solid_density, specific_body_force,
- reference_temperature};
+ intrinsic_permeability, fluid_viscosity,
+ fluid_density, biot_coefficient,
+ porosity, solid_density,
+ specific_body_force, fluid_compressibility,
+ reference_temperature, specific_gas_constant,
+ fluid_type};
SecondaryVariableCollection secondary_variables;
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
index 457acd7770b..bc6fb267bae 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
@@ -164,7 +164,11 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
pressure_size, displacement_size>::Zero(pressure_size,
displacement_size);
+ MaterialLib::Solids::MechanicsBase<DisplacementDim> const& solid_material =
+ *_process_data.solid_materials[0];
+
double const& dt = _process_data.dt;
+ double const T_ref = _process_data.reference_temperature;
ParameterLib::SpatialPosition x_position;
x_position.setElementID(_element.getID());
@@ -197,13 +201,19 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
auto& eps = _ip_data[ip].eps;
auto const& sigma_eff = _ip_data[ip].sigma_eff;
- double const S = _process_data.specific_storage(t, x_position)[0];
double const K_over_mu =
_process_data.intrinsic_permeability(t, x_position)[0] /
_process_data.fluid_viscosity(t, x_position)[0];
auto const alpha = _process_data.biot_coefficient(t, x_position)[0];
+ auto const K_S = solid_material.getBulkModulus(t, x_position);
auto const rho_sr = _process_data.solid_density(t, x_position)[0];
- auto const rho_fr = _process_data.fluid_density(t, x_position)[0];
+ // TODO (FZill) get fluid properties from GPML
+ double const p_fr =
+ (_process_data.fluid_type == FluidType::Fluid_Type::IDEAL_GAS)
+ ? N_p.dot(p)
+ : std::numeric_limits<double>::quiet_NaN();
+ double const rho_fr = _process_data.getFluidDensity(t, x_position, p_fr);
+ double const beta_p = _process_data.getFluidCompressibility(p_fr);
auto const porosity = _process_data.porosity(t, x_position)[0];
auto const& b = _process_data.specific_body_force;
auto const& identity2 = MathLib::KelvinVector::Invariants<
@@ -215,8 +225,8 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
//
eps.noalias() = B * u;
- auto C = _ip_data[ip].updateConstitutiveRelation(
- t, x_position, dt, u, _process_data.reference_temperature);
+ auto C =
+ _ip_data[ip].updateConstitutiveRelation(t, x_position, dt, u, T_ref);
local_Jac
.template block<displacement_size, displacement_size>(
@@ -239,7 +249,9 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
laplace_p.noalias() +=
rho_fr * dNdx_p.transpose() * K_over_mu * dNdx_p * w;
- storage_p.noalias() += rho_fr * N_p.transpose() * S * N_p * w;
+ storage_p.noalias() +=
+ rho_fr * N_p.transpose() * N_p * w *
+ ((alpha - porosity) * (1.0 - alpha) / K_S + porosity * beta_p);
local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
dNdx_p.transpose() * rho_fr * rho_fr * K_over_mu * b * w;
@@ -373,6 +385,9 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
pressure_size);
+ MaterialLib::Solids::MechanicsBase<DisplacementDim> const& solid_material =
+ *_process_data.solid_materials[0];
+
double const& dt = _process_data.dt;
ParameterLib::SpatialPosition x_position;
@@ -390,16 +405,25 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
auto const& N_p = _ip_data[ip].N_p;
auto const& dNdx_p = _ip_data[ip].dNdx_p;
- double const S = _process_data.specific_storage(t, x_position)[0];
double const K_over_mu =
_process_data.intrinsic_permeability(t, x_position)[0] /
_process_data.fluid_viscosity(t, x_position)[0];
auto const alpha_b = _process_data.biot_coefficient(t, x_position)[0];
- auto const rho_fr = _process_data.fluid_density(t, x_position)[0];
+ // TODO (FZill) get fluid properties from GPML
+ double const p_fr =
+ (_process_data.fluid_type == FluidType::Fluid_Type::IDEAL_GAS)
+ ? N_p.dot(p)
+ : std::numeric_limits<double>::quiet_NaN();
+ double const rho_fr = _process_data.getFluidDensity(t, x_position, p_fr);
+ double const beta_p = _process_data.getFluidCompressibility(p_fr);
+ auto const porosity = _process_data.porosity(t, x_position)[0];
+ auto const K_S = solid_material.getBulkModulus(t, x_position);
laplace.noalias() += dNdx_p.transpose() * K_over_mu * dNdx_p * w;
- mass.noalias() += N_p.transpose() * S * N_p * w;
+ mass.noalias() +=
+ N_p.transpose() * N_p * w *
+ ((alpha_b - porosity) * (1.0 - alpha_b) / K_S + porosity * beta_p);
auto const& b = _process_data.specific_body_force;
local_rhs.noalias() += dNdx_p.transpose() * rho_fr * K_over_mu * b * w;
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsProcessData.h b/ProcessLib/HydroMechanics/HydroMechanicsProcessData.h
index 7299cb661e6..9bb3ff3ff0c 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsProcessData.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsProcessData.h
@@ -10,6 +10,7 @@
#pragma once
#include "ParameterLib/Parameter.h"
+#include "MaterialLib/Fluid/FluidType/FluidType.h"
#include <memory>
#include <utility>
@@ -42,9 +43,6 @@ struct HydroMechanicsProcessData
/// Permeability of the solid. A scalar quantity,
/// ParameterLib::Parameter<double>.
ParameterLib::Parameter<double> const& intrinsic_permeability;
- /// Volumetric average specific storage of the solid and fluid phases.
- /// A scalar quantity, ParameterLib::Parameter<double>.
- ParameterLib::Parameter<double> const& specific_storage;
/// Fluid's viscosity. A scalar quantity, ParameterLib::Parameter<double>.
ParameterLib::Parameter<double> const& fluid_viscosity;
/// Fluid's density. A scalar quantity, ParameterLib::Parameter<double>.
@@ -60,9 +58,61 @@ struct HydroMechanicsProcessData
/// It is usually used to apply gravitational forces.
/// A vector of displacement dimension's length.
Eigen::Matrix<double, DisplacementDim, 1> const specific_body_force;
+
+ /// Fluid's compressibility. A scalar quantity.
+ /// Only used for compressible_fluid fluid_type
+ double const fluid_compressibility =
+ std::numeric_limits<double>::quiet_NaN();
+
+ /// Reference Temperature. A scalar quantity.
+ /// Only used for ideal_gas fluid_type
double const reference_temperature =
std::numeric_limits<double>::quiet_NaN();
+ /// Specific gas constant. A scalar quantity.
+ /// Only used for ideal_gas fluid_type
+ double const specific_gas_constant =
+ std::numeric_limits<double>::quiet_NaN();
+
+ /// Fluid type. Enumerator with possible values:
+ /// incompressible_fluid, compressible_fluid, ideal_gas
+ FluidType::Fluid_Type const fluid_type;
+
+ /// will be removed after linking with MPL
+ double getFluidDensity(
+ double const& t, ParameterLib::SpatialPosition const& x_position,
+ double const& p_fr)
+ {
+ if (fluid_type == FluidType::Fluid_Type::INCOMPRESSIBLE_FLUID ||
+ fluid_type == FluidType::Fluid_Type::COMPRESSIBLE_FLUID)
+ {
+ return fluid_density(t, x_position)[0];
+ }
+ if (fluid_type == FluidType::Fluid_Type::IDEAL_GAS)
+ {
+ return p_fr / (specific_gas_constant * reference_temperature);
+ }
+ OGS_FATAL("unknown fluid type %d", static_cast<int> (fluid_type));
+ }
+
+ /// will be removed after linking with MPL
+ double getFluidCompressibility(double const& p_fr)
+ {
+ if (fluid_type == FluidType::Fluid_Type::INCOMPRESSIBLE_FLUID)
+ {
+ return 0.0;
+ }
+ if (fluid_type == FluidType::Fluid_Type::COMPRESSIBLE_FLUID)
+ {
+ return fluid_compressibility;
+ }
+ if (fluid_type == FluidType::Fluid_Type::IDEAL_GAS)
+ {
+ return 1.0 / p_fr;
+ }
+ OGS_FATAL("unknown fluid type %d", static_cast<int> (fluid_type));
+ }
+
MeshLib::PropertyVector<double>* pressure_interpolated = nullptr;
double dt = std::numeric_limits<double>::quiet_NaN();
double t = std::numeric_limits<double>::quiet_NaN();
--
GitLab