diff --git a/ProcessLib/LiquidFlow/CreateLiquidFlowMaterialProperties.cpp b/ProcessLib/LiquidFlow/CreateLiquidFlowMaterialProperties.cpp
index 7766fbb107c04226a8941342db85504df0a52287..487d0c55796fbf1c1c7810b804f626e06a595331 100644
--- a/ProcessLib/LiquidFlow/CreateLiquidFlowMaterialProperties.cpp
+++ b/ProcessLib/LiquidFlow/CreateLiquidFlowMaterialProperties.cpp
@@ -105,12 +105,17 @@ createLiquidFlowMaterialProperties(
"thermal_expansion", parameters, 1);
DBUG("Use \'%s\' as solid thermal expansion.",
solid_thermal_expansion.name.c_str());
+ auto& biot_constant = findParameter<double>(
+ *solid_config,
+ //! \ogs_file_param{prj__processes__process__LIQUID_FLOW__material_property__solid__biot_constant}
+ "biot_constant", parameters, 1);
return std::unique_ptr<LiquidFlowMaterialProperties>(
new LiquidFlowMaterialProperties(
std::move(fluid_properties),
std::move(intrinsic_permeability_models),
std::move(porosity_models), std::move(storage_models),
- has_material_ids, material_ids, solid_thermal_expansion));
+ has_material_ids, material_ids, solid_thermal_expansion,
+ biot_constant));
}
else
{
@@ -121,7 +126,8 @@ createLiquidFlowMaterialProperties(
std::move(fluid_properties),
std::move(intrinsic_permeability_models),
std::move(porosity_models), std::move(storage_models),
- has_material_ids, material_ids, void_parameter));
+ has_material_ids, material_ids, void_parameter,
+ void_parameter));
}
}
diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
index 1cedbbb231e7b356058621f45f396a1728504ed2..c28b5e44fdb05404ce5d4f5d9212a6ced04f7204 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
@@ -60,20 +60,45 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
{
SpatialPosition pos;
pos.setElementID(_element.getID());
- _material_properties->setMaterialID(pos);
+ const int material_id = _material_properties.getMaterialID(pos);
- const Eigen::MatrixXd& perm =
- _material_properties->getPermeability(t, pos, _element.getDimension());
+ const Eigen::MatrixXd& perm = _material_properties.getPermeability(
+ material_id, t, pos, _element.getDimension());
// Note: For Inclined 1D in 2D/3D or 2D element in 3D, the first item in
// the assert must be changed to perm.rows() == _element->getDimension()
assert(perm.rows() == GlobalDim || perm.rows() == 1);
- if (perm.size() == 1) // isotropic or 1D problem.
- local_assemble<IsotropicCalculator>(
- t, local_x, local_M_data, local_K_data, local_b_data, pos, perm);
- else
- local_assemble<AnisotropicCalculator>(
- t, local_x, local_M_data, local_K_data, local_b_data, pos, perm);
+ const double dt = coupled_term.dt;
+ auto it = coupled_term.coupled_processes.begin();
+ while (it != coupled_term.coupled_processes.end())
+ {
+ switch (it->first)
+ {
+ case ProcessLib::ProcessType::HeatConductionProcess:
+ {
+ const auto local_T0 = coupled_term.local_coupled_xs0.at(
+ ProcessLib::ProcessType::HeatConductionProcess);
+ const auto local_T1 = coupled_term.local_coupled_xs.at(
+ ProcessLib::ProcessType::HeatConductionProcess);
+
+ if (perm.size() == 1) // isotropic or 1D problem.
+ local_assembleCoupledWithHeatTransport<IsotropicCalculator>(
+ material_id, t, dt, local_x, local_T0, local_T1,
+ local_M_data, local_K_data, local_b_data, pos, perm);
+ else
+ local_assembleCoupledWithHeatTransport<
+ AnisotropicCalculator>(
+ material_id, t, dt, local_x, local_T0, local_T1,
+ local_M_data, local_K_data, local_b_data, pos, perm);
+ }
+ break;
+ default:
+ OGS_FATAL(
+ "This coupled process is not presented for "
+ "LiquidFlow process");
+ }
+ it++;
+ }
}
template <typename ShapeFunction, typename IntegrationMethod,
@@ -142,8 +167,9 @@ template <typename ShapeFunction, typename IntegrationMethod,
template <typename LaplacianGravityVelocityCalculator>
void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
local_assembleCoupledWithHeatTransport(
- double const t, std::vector<double> const& local_x,
- std::vector<double> const& local_T, std::vector<double>& local_M_data,
+ const int material_id, double const t, double const dt,
+ std::vector<double> const& local_x, std::vector<double> const& local_T0,
+ std::vector<double> const& local_T1, std::vector<double>& local_M_data,
std::vector<double>& local_K_data, std::vector<double>& local_b_data,
SpatialPosition const& pos, Eigen::MatrixXd const& perm)
{
@@ -171,28 +197,44 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
double p = 0.;
NumLib::shapeFunctionInterpolate(local_x, sm.N, p);
+ double T0 = 0.;
+ NumLib::shapeFunctionInterpolate(local_T0, sm.N, T0);
double T = 0.;
- NumLib::shapeFunctionInterpolate(local_T, sm.N, T);
+ NumLib::shapeFunctionInterpolate(local_T1, sm.N, T);
const double integration_factor =
sm.integralMeasure * sm.detJ * wp.getWeight();
// Assemble mass matrix, M
- local_M.noalias() +=
- _material_properties->getMassCoefficient(t, pos, porosity_variable,
- storage_variable, p, T) *
- sm.N.transpose() * sm.N * integration_factor;
+ local_M.noalias() += _material_properties.getMassCoefficient(
+ material_id, t, pos, porosity_variable,
+ storage_variable, p, T) *
+ sm.N.transpose() * sm.N * integration_factor;
// Compute density:
- const double rho_g = _material_properties->getLiquidDensity(p, T) *
- _gravitational_acceleration;
+ const double rho = _material_properties.getLiquidDensity(p, T);
+ const double rho_g = rho * _gravitational_acceleration;
+
// Compute viscosity:
- const double mu = _material_properties->getViscosity(p, T);
+ const double mu = _material_properties.getViscosity(p, T);
// Assemble Laplacian, K, and RHS by the gravitational term
LaplacianGravityVelocityCalculator::calculateLaplacianAndGravityTerm(
local_K, local_b, sm, perm, integration_factor, mu, rho_g,
_gravitational_axis_id);
+
+ // Add the thermal expansion term
+ auto const solid_thermal_expansion =
+ _material_properties.getSolidThermalExpansion(t, pos);
+ auto const biot_constant =
+ _material_properties.getBiotConstant(t, pos);
+ auto const porosity = _material_properties.getPorosity(
+ material_id, porosity_variable, T);
+ const double eff_thermal_expansion =
+ 3.0 * (biot_constant - porosity) * solid_thermal_expansion +
+ porosity * _material_properties.getdLiquidDensity_dT(p, T) / rho;
+ local_b.noalias() +=
+ eff_thermal_expansion * (T - T0) * integration_factor * sm.N / dt;
}
}
@@ -258,7 +300,6 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
}
}
-
template <typename ShapeFunction, typename IntegrationMethod,
unsigned GlobalDim>
template <typename LaplacianGravityVelocityCalculator>
@@ -287,10 +328,10 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
double T = 0.;
NumLib::shapeFunctionInterpolate(local_T, sm.N, T);
- const double rho_g = _material_properties->getLiquidDensity(p, T) *
+ const double rho_g = _material_properties.getLiquidDensity(p, T) *
_gravitational_acceleration;
// Compute viscosity:
- const double mu = _material_properties->getViscosity(p, T);
+ const double mu = _material_properties.getViscosity(p, T);
LaplacianGravityVelocityCalculator::calculateVelocity(
_darcy_velocities, local_p_vec, sm, perm, ip, mu, rho_g,
diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
index 884fc8d2b6afba9863d04ed441a8e244a83bb32e..c8a61de086c371df3423d31cf1c92fc6bebe5e37 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
@@ -190,14 +190,11 @@ private:
template <typename LaplacianGravityVelocityCalculator>
void local_assembleCoupledWithHeatTransport(
- double const t,
- std::vector<double> const& local_x,
- std::vector<double> const& local_T,
- std::vector<double>& local_M_data,
- std::vector<double>& local_K_data,
- std::vector<double>& local_b_data,
- SpatialPosition const& pos,
- Eigen::MatrixXd const& perm);
+ const int material_id, double const t, double const dt,
+ std::vector<double> const& local_x, std::vector<double> const& local_T0,
+ std::vector<double> const& local_T1, std::vector<double>& local_M_data,
+ std::vector<double>& local_K_data, std::vector<double>& local_b_data,
+ SpatialPosition const& pos, Eigen::MatrixXd const& perm);
template <typename LaplacianGravityVelocityCalculator>
void computeSecondaryVariableLocal(double const /*t*/,
diff --git a/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.cpp b/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.cpp
index 295423d505607bd7ade397e0425ae0221608b99d..2f99c1360cff4e955f7a0a1470d7859660a92501 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.cpp
+++ b/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.cpp
@@ -55,6 +55,17 @@ double LiquidFlowMaterialProperties::getLiquidDensity(const double p,
MaterialLib::Fluid::FluidPropertyType::Density, vars);
}
+double LiquidFlowMaterialProperties::getdLiquidDensity_dT(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 _fluid_properties->getdValue(
+ MaterialLib::Fluid::FluidPropertyType::Density, vars,
+ MaterialLib::Fluid::PropertyVariableType::T);
+}
+
double LiquidFlowMaterialProperties::getViscosity(const double p,
const double T) const
{
@@ -65,6 +76,26 @@ double LiquidFlowMaterialProperties::getViscosity(const double p,
MaterialLib::Fluid::FluidPropertyType::Viscosity, vars);
}
+double LiquidFlowMaterialProperties::getHeatCapacity(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 _fluid_properties->getValue(
+ MaterialLib::Fluid::FluidPropertyType::HeatCapacity, vars);
+}
+
+double LiquidFlowMaterialProperties::getThermalConductivity(
+ 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 _fluid_properties->getValue(
+ MaterialLib::Fluid::FluidPropertyType::ThermalConductivity, vars);
+}
+
double LiquidFlowMaterialProperties::getMassCoefficient(
const int material_id, const double /*t*/, const SpatialPosition& /*pos*/,
const double p, const double T, const double porosity_variable,
@@ -94,5 +125,17 @@ Eigen::MatrixXd const& LiquidFlowMaterialProperties::getPermeability(
return _intrinsic_permeability_models[material_id];
}
+double LiquidFlowMaterialProperties::getSolidThermalExpansion(
+ const double t, const SpatialPosition& pos) const
+{
+ return _solid_thermal_expansion(t, pos)[0];
+}
+
+double LiquidFlowMaterialProperties::getBiotConstant(
+ const double t, const SpatialPosition& pos) const
+{
+ return _biot_constant(t, pos)[0];
+}
+
} // end of namespace
} // end of namespace
diff --git a/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.h b/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.h
index 2078f1608ddb8f97e4cd9757a42ef57d839eca9e..90f5557e835eeb23e2810b8922d24748c49e62b6 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.h
@@ -68,7 +68,8 @@ public:
storage_models,
bool const has_material_ids,
MeshLib::PropertyVector<int> const& material_ids,
- Parameter<double> const& solid_thermal_expansion)
+ Parameter<double> const& solid_thermal_expansion,
+ Parameter<double> const& biot_constant)
: _has_material_ids(has_material_ids),
_material_ids(material_ids),
_fluid_properties(std::move(fluid_properties)),
@@ -76,7 +77,8 @@ public:
std::move(intrinsic_permeability_models)),
_porosity_models(std::move(porosity_models)),
_storage_models(std::move(storage_models)),
- _solid_thermal_expansion(solid_thermal_expansion)
+ _solid_thermal_expansion(solid_thermal_expansion),
+ _biot_constant(biot_constant)
{
}
@@ -111,8 +113,25 @@ public:
double getLiquidDensity(const double p, const double T) const;
+ double getdLiquidDensity_dT(const double p, const double T) const;
+
double getViscosity(const double p, const double T) const;
+ double getHeatCapacity(const double p, const double T) const;
+
+ double getThermalConductivity(const double p, const double T) const;
+
+ double getPorosity(const int material_id, const double porosity_variable,
+ const double T) const
+ {
+ return _porosity_models[material_id]->getValue(porosity_variable, T);
+ }
+
+ double getSolidThermalExpansion(const double t,
+ const SpatialPosition& pos) const;
+
+ double getBiotConstant(const double t, const SpatialPosition& pos) const;
+
private:
/// A flag to indicate whether the reference member, _material_ids,
/// is not assigned.
@@ -132,6 +151,7 @@ private:
_storage_models;
Parameter<double> const& _solid_thermal_expansion;
+ Parameter<double> const& _biot_constant;
// Note: For the statistical data of porous media, they could be read from
// vtu files directly. This can be done by using property vectors directly.
// Such property vectors will be added here if they are needed.