diff --git a/ProcessLib/HT/HTFEM.h b/ProcessLib/HT/HTFEM.h
index c3ae25856f5268e050f1a390f9df624437d3a662..da24a85a76705cee7c867691e7e4806e9880ec8d 100644
--- a/ProcessLib/HT/HTFEM.h
+++ b/ProcessLib/HT/HTFEM.h
@@ -162,6 +162,56 @@ protected:
Eigen::aligned_allocator<
IntegrationPointData<NodalRowVectorType, GlobalDimNodalMatrixType>>>
_ip_data;
+
+ double getHeatEnergyCoefficient(const double t, const SpatialPosition& pos,
+ const double porosity,
+ const double fluid_density,
+ const double specific_heat_capacity_fluid)
+ {
+ auto const& material_properties = this->_material_properties;
+ auto const specific_heat_capacity_solid =
+ material_properties.specific_heat_capacity_solid(t, pos)[0];
+
+ auto const solid_density = material_properties.density_solid(t, pos)[0];
+
+ return solid_density * specific_heat_capacity_solid * (1 - porosity) +
+ fluid_density * specific_heat_capacity_fluid * porosity;
+ }
+
+ GlobalDimMatrixType getThermalConductivityDispersivity(
+ const double t, const SpatialPosition& pos, const double porosity,
+ const double fluid_density, const double specific_heat_capacity_fluid,
+ const GlobalDimVectorType& velocity, const GlobalDimMatrixType& I)
+ {
+ auto const& material_properties = this->_material_properties;
+
+ auto const thermal_conductivity_solid =
+ material_properties.thermal_conductivity_solid(t, pos)[0];
+ auto const thermal_conductivity_fluid =
+ material_properties.thermal_conductivity_fluid(t, pos)[0];
+ double const thermal_conductivity =
+ thermal_conductivity_solid * (1 - porosity) +
+ thermal_conductivity_fluid * porosity;
+
+ auto const thermal_dispersivity_longitudinal =
+ material_properties.thermal_dispersivity_longitudinal(t, pos)[0];
+ auto const thermal_dispersivity_transversal =
+ material_properties.thermal_dispersivity_transversal(t, pos)[0];
+
+ if (thermal_dispersivity_longitudinal == 0.0 &&
+ thermal_dispersivity_transversal == 0.0)
+ return thermal_conductivity * I;
+
+ double const velocity_magnitude = velocity.norm();
+ GlobalDimMatrixType const thermal_dispersivity =
+ fluid_density * specific_heat_capacity_fluid *
+ (thermal_dispersivity_transversal * velocity_magnitude * I +
+ (thermal_dispersivity_longitudinal -
+ thermal_dispersivity_transversal) /
+ velocity_magnitude * velocity * velocity.transpose());
+
+ return thermal_conductivity * I + thermal_dispersivity;
+ }
};
} // namespace HT
diff --git a/ProcessLib/HT/HTMaterialProperties.h b/ProcessLib/HT/HTMaterialProperties.h
index 5cea0c469433d4b5505835552a12cb2a99e35d3c..57d96f236af26e1e63195d459c3ab66edc29b3e6 100644
--- a/ProcessLib/HT/HTMaterialProperties.h
+++ b/ProcessLib/HT/HTMaterialProperties.h
@@ -56,15 +56,6 @@ struct HTMaterialProperties
{
}
- //! Copies are forbidden.
- HTMaterialProperties(HTMaterialProperties const&) = delete;
-
- //! Assignments are not needed.
- void operator=(HTMaterialProperties const&) = delete;
-
- //! Assignments are not needed.
- void operator=(HTMaterialProperties&&) = delete;
-
MaterialLib::PorousMedium::PorousMediaProperties porous_media_properties;
Parameter<double> const& density_solid;
Parameter<double> const& fluid_reference_density;
diff --git a/ProcessLib/HT/MonolithicHTFEM.h b/ProcessLib/HT/MonolithicHTFEM.h
index 56e8bb55c0fdfb7a2e6d996a8b6e6af521199c05..6aed470778d09b07eda59d2a7017182a8e26c89c 100644
--- a/ProcessLib/HT/MonolithicHTFEM.h
+++ b/ProcessLib/HT/MonolithicHTFEM.h
@@ -112,21 +112,12 @@ public:
{
pos.setIntegrationPoint(ip);
- auto const fluid_reference_density =
- process_data.fluid_reference_density(t, pos)[0];
-
- auto const density_solid = process_data.density_solid(t, pos)[0];
// \todo the argument to getValue() has to be changed for non
// constant storage model
auto const specific_storage =
process_data.porous_media_properties.getSpecificStorage(t, pos)
.getValue(0.0);
- auto const thermal_conductivity_solid =
- process_data.thermal_conductivity_solid(t, pos)[0];
- auto const thermal_conductivity_fluid =
- process_data.thermal_conductivity_fluid(t, pos)[0];
-
auto const& ip_data = this->_ip_data[ip];
auto const& N = ip_data.N;
auto const& dNdx = ip_data.dNdx;
@@ -146,13 +137,6 @@ public:
process_data.porous_media_properties.getIntrinsicPermeability(
t, pos).getValue(t, pos, 0.0, T_int_pt);
- double const thermal_conductivity =
- thermal_conductivity_solid * (1 - porosity) +
- thermal_conductivity_fluid * porosity;
-
- auto const specific_heat_capacity_solid =
- process_data.specific_heat_capacity_solid(t, pos)[0];
-
vars[static_cast<int>(
MaterialLib::Fluid::PropertyVariableType::T)] = T_int_pt;
vars[static_cast<int>(
@@ -161,13 +145,8 @@ public:
process_data.fluid_properties->getValue(
MaterialLib::Fluid::FluidPropertyType::HeatCapacity, vars);
- auto const thermal_dispersivity_longitudinal =
- process_data.thermal_dispersivity_longitudinal(t, pos)[0];
- auto const thermal_dispersivity_transversal =
- process_data.thermal_dispersivity_transversal(t, pos)[0];
-
// Use the fluid density model to compute the density
- auto const density = process_data.fluid_properties->getValue(
+ auto const fluid_density = process_data.fluid_properties->getValue(
MaterialLib::Fluid::FluidPropertyType::Density, vars);
// Use the viscosity model to compute the viscosity
@@ -177,37 +156,29 @@ public:
GlobalDimVectorType const velocity =
process_data.has_gravity
- ? GlobalDimVectorType(-K_over_mu *
- (dNdx * p_nodal_values - density * b))
+ ? GlobalDimVectorType(-K_over_mu * (dNdx * p_nodal_values -
+ fluid_density * b))
: GlobalDimVectorType(-K_over_mu * dNdx * p_nodal_values);
- double const velocity_magnitude = velocity.norm();
- GlobalDimMatrixType const thermal_dispersivity =
- fluid_reference_density * specific_heat_capacity_fluid *
- (thermal_dispersivity_transversal * velocity_magnitude * I +
- (thermal_dispersivity_longitudinal -
- thermal_dispersivity_transversal) /
- velocity_magnitude * velocity * velocity.transpose());
-
- GlobalDimMatrixType const hydrodynamic_thermodispersion =
- thermal_conductivity * I + thermal_dispersivity;
-
- double const heat_capacity =
- density_solid * specific_heat_capacity_solid * (1 - porosity) +
- fluid_reference_density * specific_heat_capacity_fluid *
- porosity;
-
// matrix assembly
+ GlobalDimMatrixType const thermal_conductivity_dispersivity =
+ this->getThermalConductivityDispersivity(
+ t, pos, porosity, fluid_density,
+ specific_heat_capacity_fluid, velocity, I);
Ktt.noalias() +=
- (dNdx.transpose() * hydrodynamic_thermodispersion * dNdx +
- N.transpose() * velocity.transpose() * dNdx *
- fluid_reference_density * specific_heat_capacity_fluid) *
+ (dNdx.transpose() * thermal_conductivity_dispersivity * dNdx +
+ N.transpose() * velocity.transpose() * dNdx * fluid_density *
+ specific_heat_capacity_fluid) *
w;
Kpp.noalias() += w * dNdx.transpose() * K_over_mu * dNdx;
- Mtt.noalias() += w * N.transpose() * heat_capacity * N;
+ Mtt.noalias() +=
+ w *
+ this->getHeatEnergyCoefficient(t, pos, porosity, fluid_density,
+ specific_heat_capacity_fluid) *
+ N.transpose() * N;
Mpp.noalias() += w * N.transpose() * specific_storage * N;
if (process_data.has_gravity)
- Bp += w * density * dNdx.transpose() * K_over_mu * b;
+ Bp += w * fluid_density * dNdx.transpose() * K_over_mu * b;
/* with Oberbeck-Boussing assumption density difference only exists
* in buoyancy effects */
}
diff --git a/ProcessLib/HT/StaggeredHTFEM-impl.h b/ProcessLib/HT/StaggeredHTFEM-impl.h
index ec6727ad42cd895404c99e4900bcf6cad745aa6c..5a2a916180d19d02a4b2f664a4327728c753f448 100644
--- a/ProcessLib/HT/StaggeredHTFEM-impl.h
+++ b/ProcessLib/HT/StaggeredHTFEM-impl.h
@@ -225,44 +225,21 @@ void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
T1_at_xi;
vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::p)] =
p_at_xi;
-
- // Assemble mass matrix
- auto const specific_heat_capacity_solid =
- material_properties.specific_heat_capacity_solid(t, pos)[0];
- auto const specific_heat_capacity_fluid =
- material_properties.fluid_properties->getValue(
- MaterialLib::Fluid::FluidPropertyType::HeatCapacity, vars);
- auto const fluid_reference_density =
- material_properties.fluid_reference_density(t, pos)[0];
-
- auto const solid_density = material_properties.density_solid(t, pos)[0];
-
- // Use the fluid density model to compute the density
auto const fluid_density =
material_properties.fluid_properties->getValue(
MaterialLib::Fluid::FluidPropertyType::Density, vars);
- double const heat_capacity =
- solid_density * specific_heat_capacity_solid * (1 - porosity) +
- fluid_density * specific_heat_capacity_fluid * porosity;
+ auto const specific_heat_capacity_fluid =
+ material_properties.fluid_properties->getValue(
+ MaterialLib::Fluid::FluidPropertyType::HeatCapacity, vars);
- local_M.noalias() += w * heat_capacity * N.transpose() * N;
+ // Assemble mass matrix
+ local_M.noalias() +=
+ w *
+ this->getHeatEnergyCoefficient(t, pos, porosity, fluid_density,
+ specific_heat_capacity_fluid) *
+ N.transpose() * N;
// Assemble Laplace matrix
-
- auto const thermal_conductivity_solid =
- material_properties.thermal_conductivity_solid(t, pos)[0];
- auto const thermal_conductivity_fluid =
- material_properties.thermal_conductivity_fluid(t, pos)[0];
- double const thermal_conductivity =
- thermal_conductivity_solid * (1 - porosity) +
- thermal_conductivity_fluid * porosity;
-
- auto const thermal_dispersivity_longitudinal =
- material_properties.thermal_dispersivity_longitudinal(t, pos)[0];
- auto const thermal_dispersivity_transversal =
- material_properties.thermal_dispersivity_transversal(t, pos)[0];
-
- // Use the viscosity model to compute the viscosity
auto const viscosity = material_properties.fluid_properties->getValue(
MaterialLib::Fluid::FluidPropertyType::Viscosity, vars);
@@ -278,19 +255,13 @@ void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
fluid_density * b))
: GlobalDimVectorType(-K_over_mu * dNdx * local_p_Eigen_type);
- double const velocity_magnitude = velocity.norm();
- GlobalDimMatrixType const thermal_dispersivity =
- fluid_reference_density * specific_heat_capacity_fluid *
- (thermal_dispersivity_transversal * velocity_magnitude * I +
- (thermal_dispersivity_longitudinal -
- thermal_dispersivity_transversal) /
- velocity_magnitude * velocity * velocity.transpose());
-
- GlobalDimMatrixType const hydrodynamic_thermodispersion =
- thermal_conductivity * I + thermal_dispersivity;
+ GlobalDimMatrixType const thermal_conductivity_dispersivity =
+ this->getThermalConductivityDispersivity(
+ t, pos, porosity, fluid_density, specific_heat_capacity_fluid,
+ velocity, I);
local_K.noalias() +=
- w * (dNdx.transpose() * hydrodynamic_thermodispersion * dNdx +
+ w * (dNdx.transpose() * thermal_conductivity_dispersivity * dNdx +
N.transpose() * velocity.transpose() * dNdx * fluid_density *
specific_heat_capacity_fluid);
}