diff --git a/ProcessLib/HT/HTFEM.h b/ProcessLib/HT/HTFEM.h
index 561dd39acef7ed6763072c69662189d96621cb3f..c3ae25856f5268e050f1a390f9df624437d3a662 100644
--- a/ProcessLib/HT/HTFEM.h
+++ b/ProcessLib/HT/HTFEM.h
@@ -47,10 +47,10 @@ public:
std::size_t const local_matrix_size,
bool is_axially_symmetric,
unsigned const integration_order,
- HTMaterialProperties const& process_data,
+ HTMaterialProperties const& material_properties,
const unsigned dof_per_node)
: _element(element),
- _process_data(process_data),
+ _material_properties(material_properties),
_integration_method(integration_order)
{
// This assertion is valid only if all nodal d.o.f. use the same shape
@@ -129,20 +129,21 @@ public:
MaterialLib::Fluid::PropertyVariableType::p)] = p_int_pt;
auto const K =
- _process_data.porous_media_properties.getIntrinsicPermeability(
+ _material_properties.porous_media_properties.getIntrinsicPermeability(
t, pos).getValue(t, pos, 0.0, T_int_pt);
- auto const mu = _process_data.fluid_properties->getValue(
+ auto const mu = _material_properties.fluid_properties->getValue(
MaterialLib::Fluid::FluidPropertyType::Viscosity, vars);
GlobalDimMatrixType const K_over_mu = K / mu;
cache_mat.col(ip).noalias() = -K_over_mu * dNdx * p_nodal_values;
- if (_process_data.has_gravity)
+ if (_material_properties.has_gravity)
{
- auto const rho_w = _process_data.fluid_properties->getValue(
- MaterialLib::Fluid::FluidPropertyType::Density, vars);
- auto const b = _process_data.specific_body_force;
+ auto const rho_w =
+ _material_properties.fluid_properties->getValue(
+ MaterialLib::Fluid::FluidPropertyType::Density, vars);
+ auto const b = _material_properties.specific_body_force;
// here it is assumed that the vector b is directed 'downwards'
cache_mat.col(ip).noalias() += K_over_mu * rho_w * b;
}
@@ -153,7 +154,7 @@ public:
protected:
MeshLib::Element const& _element;
- HTMaterialProperties const& _process_data;
+ HTMaterialProperties const& _material_properties;
IntegrationMethod const _integration_method;
std::vector<
diff --git a/ProcessLib/HT/HTMaterialProperties.h b/ProcessLib/HT/HTMaterialProperties.h
index 9f71a1ddb15ece2ffa1611396b44305f6e25bfc0..5cea0c469433d4b5505835552a12cb2a99e35d3c 100644
--- a/ProcessLib/HT/HTMaterialProperties.h
+++ b/ProcessLib/HT/HTMaterialProperties.h
@@ -50,29 +50,12 @@ struct HTMaterialProperties
thermal_conductivity_solid(thermal_conductivity_solid_),
thermal_conductivity_fluid(thermal_conductivity_fluid_),
solid_thermal_expansion(solid_thermal_expansion_),
- biot_constant(biot_constant),
+ biot_constant(biot_constant_),
specific_body_force(std::move(specific_body_force_)),
has_gravity(has_gravity_)
{
}
- HTMaterialProperties(HTMaterialProperties&& other)
- : porous_media_properties(std::move(other.porous_media_properties)),
- density_solid(other.density_solid),
- fluid_reference_density(other.fluid_reference_density),
- fluid_properties(other.fluid_properties.release()),
- specific_heat_capacity_solid(other.specific_heat_capacity_solid),
- thermal_dispersivity_longitudinal(
- other.thermal_dispersivity_longitudinal),
- thermal_dispersivity_transversal(
- other.thermal_dispersivity_transversal),
- thermal_conductivity_solid(other.thermal_conductivity_solid),
- thermal_conductivity_fluid(other.thermal_conductivity_fluid),
- specific_body_force(other.specific_body_force),
- has_gravity(other.has_gravity)
- {
- }
-
//! Copies are forbidden.
HTMaterialProperties(HTMaterialProperties const&) = delete;
diff --git a/ProcessLib/HT/HTProcess.cpp b/ProcessLib/HT/HTProcess.cpp
index 3d993042da7effe60df5cf170a569d51148fab14..c17570c121a9d2a24a3068b6891a8eb1abe473b6 100644
--- a/ProcessLib/HT/HTProcess.cpp
+++ b/ProcessLib/HT/HTProcess.cpp
@@ -50,14 +50,16 @@ void HTProcess::initializeConcreteProcess(
ProcessLib::createLocalAssemblers<MonolithicHTFEM>(
mesh.getDimension(), mesh.getElements(), dof_table,
pv.getShapeFunctionOrder(), _local_assemblers,
- mesh.isAxiallySymmetric(), integration_order, *_material_properties);
+ mesh.isAxiallySymmetric(), integration_order,
+ *_material_properties);
}
else
{
ProcessLib::createLocalAssemblers<StaggeredHTFEM>(
mesh.getDimension(), mesh.getElements(), dof_table,
pv.getShapeFunctionOrder(), _local_assemblers,
- mesh.isAxiallySymmetric(), integration_order, *_material_properties);
+ mesh.isAxiallySymmetric(), integration_order,
+ *_material_properties);
}
_secondary_variables.addSecondaryVariable(
diff --git a/ProcessLib/HT/HTProcess.h b/ProcessLib/HT/HTProcess.h
index ce773e5807b0dfc8a189d56c6e9746152ee0e419..61b399d572a117b0af4f2eeddb62c339aa0b7cab 100644
--- a/ProcessLib/HT/HTProcess.h
+++ b/ProcessLib/HT/HTProcess.h
@@ -9,6 +9,8 @@
#pragma once
+#include <array>
+
#include "NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h"
#include "ProcessLib/Process.h"
@@ -87,10 +89,9 @@ private:
const std::unique_ptr<HTMaterialProperties> _material_properties;
std::vector<std::unique_ptr<HTLocalAssemblerInterface>> _local_assemblers;
-
+
/// Solutions of the previous time step
std::array<std::unique_ptr<GlobalVector>, 2> _xs_previous_timestep;
-
};
} // namespace HT
diff --git a/ProcessLib/HT/MonolithicHTFEM.h b/ProcessLib/HT/MonolithicHTFEM.h
index 04d167672d2b5d0f0fcf75b647c5bb27ee9fdeef..56e8bb55c0fdfb7a2e6d996a8b6e6af521199c05 100644
--- a/ProcessLib/HT/MonolithicHTFEM.h
+++ b/ProcessLib/HT/MonolithicHTFEM.h
@@ -56,10 +56,10 @@ public:
std::size_t const local_matrix_size,
bool is_axially_symmetric,
unsigned const integration_order,
- HTMaterialProperties const& process_data)
+ HTMaterialProperties const& material_properties)
: HTFEM<ShapeFunction, IntegrationMethod, GlobalDim>(
element, local_matrix_size, is_axially_symmetric,
- integration_order, process_data, NUM_NODAL_DOF)
+ integration_order, material_properties, NUM_NODAL_DOF)
{
}
@@ -96,7 +96,7 @@ public:
auto p_nodal_values =
Eigen::Map<const NodalVectorType>(&local_x[num_nodes], num_nodes);
- auto const& process_data = this->_process_data;
+ auto const& process_data = this->_material_properties;
auto const& b = process_data.specific_body_force;
diff --git a/ProcessLib/HT/StaggeredHTFEM-impl.h b/ProcessLib/HT/StaggeredHTFEM-impl.h
index 1d6508049c4f73843c920d67ea231e8c5d0da37a..ec6727ad42cd895404c99e4900bcf6cad745aa6c 100644
--- a/ProcessLib/HT/StaggeredHTFEM-impl.h
+++ b/ProcessLib/HT/StaggeredHTFEM-impl.h
@@ -48,15 +48,15 @@ void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
LocalCoupledSolutions const& coupled_term)
{
auto const& local_p = coupled_term.local_coupled_xs[0];
- auto const& local_T = coupled_term.local_coupled_xs[1];
- auto const& local_T1 = coupled_term.local_coupled_xs0[1];
- const double dt = coupled_term.dt;
-
- auto const local_matrix_size = local_x.size();
+ auto const local_matrix_size = local_p.size();
// This assertion is valid only if all nodal d.o.f. use the same shape
// matrices.
assert(local_matrix_size == ShapeFunction::NPOINTS);
+ auto const& local_T1 = coupled_term.local_coupled_xs[1];
+ auto const& local_T0 = coupled_term.local_coupled_xs0[1];
+ const double dt = coupled_term.dt;
+
auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
local_M_data, local_matrix_size, local_matrix_size);
auto local_K = MathLib::createZeroedMatrix<LocalMatrixType>(
@@ -67,9 +67,9 @@ void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
SpatialPosition pos;
pos.setElementID(this->_element.getID());
- auto const& process_data = this->_process_data;
+ auto const& material_properties = this->_material_properties;
- auto const& b = process_data.specific_body_force;
+ auto const& b = material_properties.specific_body_force;
GlobalDimMatrixType const& I(
GlobalDimMatrixType::Identity(GlobalDim, GlobalDim));
@@ -83,74 +83,200 @@ void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
{
pos.setIntegrationPoint(ip);
- auto const fluid_reference_density =
- process_data.fluid_reference_density(t, pos)[0];
+ auto const& ip_data = this->_ip_data[ip];
+ auto const& N = ip_data.N;
+ auto const& dNdx = ip_data.dNdx;
+ auto const& w = ip_data.integration_weight;
+
+ double p_at_xi = 0.;
+ NumLib::shapeFunctionInterpolate(local_p, N, p_at_xi);
+ double T1_at_xi = 0.;
+ NumLib::shapeFunctionInterpolate(local_T1, N, T1_at_xi);
+
+ auto const porosity =
+ material_properties.porous_media_properties.getPorosity(t, pos)
+ .getValue(t, pos, 0.0, T1_at_xi);
+
+ vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)] =
+ T1_at_xi;
+ vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::p)] =
+ p_at_xi;
+
+ // Use the fluid density model to compute the density
+ auto const fluid_density =
+ material_properties.fluid_properties->getValue(
+ MaterialLib::Fluid::FluidPropertyType::Density, vars);
+ const double dfluid_density_dp =
+ material_properties.fluid_properties->getdValue(
+ MaterialLib::Fluid::FluidPropertyType::Density, vars,
+ MaterialLib::Fluid::PropertyVariableType::p);
+ const double dfluid_density_dT =
+ material_properties.fluid_properties->getdValue(
+ MaterialLib::Fluid::FluidPropertyType::Density, vars,
+ MaterialLib::Fluid::PropertyVariableType::T);
+
+ // Use the viscosity model to compute the viscosity
+ auto const viscosity = material_properties.fluid_properties->getValue(
+ MaterialLib::Fluid::FluidPropertyType::Viscosity, vars);
- 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)
+ material_properties.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 intrinsic_permeability =
+ material_properties.porous_media_properties
+ .getIntrinsicPermeability(t, pos)
+ .getValue(t, pos, 0.0, T1_at_xi);
+ GlobalDimMatrixType K_over_mu = intrinsic_permeability / viscosity;
+
+ // matrix assembly
+ local_M.noalias() += w * (porosity * dfluid_density_dp / fluid_density +
+ specific_storage) *
+ N.transpose() * N;
+
+ local_K.noalias() += w * dNdx.transpose() * K_over_mu * dNdx;
+
+ // Add the thermal expansion term
+ auto const solid_thermal_expansion =
+ material_properties.solid_thermal_expansion(t, pos)[0];
+ if (solid_thermal_expansion > 0.0)
+ {
+ double T0_at_xi = 0.;
+ NumLib::shapeFunctionInterpolate(local_T0, N, T0_at_xi);
+ auto const biot_constant =
+ material_properties.biot_constant(t, pos)[0];
+ const double eff_thermal_expansion =
+ 3.0 * (biot_constant - porosity) * solid_thermal_expansion -
+ porosity * dfluid_density_dT / fluid_density;
+ local_b.noalias() +=
+ eff_thermal_expansion * (T1_at_xi - T0_at_xi) * w * N / dt;
+ }
+
+ if (material_properties.has_gravity)
+ {
+ local_b.noalias() +=
+ w * fluid_density * dNdx.transpose() * K_over_mu * b;
+ }
+ }
+}
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
+ assembleHeatTransportEquation(double const t,
+ std::vector<double>& local_M_data,
+ std::vector<double>& local_K_data,
+ std::vector<double>& /*local_b_data*/,
+ LocalCoupledSolutions const& coupled_term)
+{
+ auto const& local_p = coupled_term.local_coupled_xs[0];
+ auto const local_matrix_size = local_p.size();
+ // This assertion is valid only if all nodal d.o.f. use the same shape
+ // matrices.
+ assert(local_matrix_size == ShapeFunction::NPOINTS);
+
+ auto local_p_Eigen_type =
+ Eigen::Map<const NodalVectorType>(&local_p[0], local_matrix_size);
+
+ auto const& local_T1 = coupled_term.local_coupled_xs[1];
+
+ auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
+ local_M_data, local_matrix_size, local_matrix_size);
+ auto local_K = MathLib::createZeroedMatrix<LocalMatrixType>(
+ local_K_data, local_matrix_size, local_matrix_size);
+
+ SpatialPosition pos;
+ pos.setElementID(this->_element.getID());
+
+ auto const& material_properties = this->_material_properties;
+
+ auto const& b = material_properties.specific_body_force;
+
+ GlobalDimMatrixType const& I(
+ GlobalDimMatrixType::Identity(GlobalDim, GlobalDim));
+
+ MaterialLib::Fluid::FluidProperty::ArrayType vars;
+
+ unsigned const n_integration_points =
+ this->_integration_method.getNumberOfPoints();
+
+ for (std::size_t ip(0); ip < n_integration_points; ip++)
+ {
+ pos.setIntegrationPoint(ip);
auto const& ip_data = this->_ip_data[ip];
auto const& N = ip_data.N;
auto const& dNdx = ip_data.dNdx;
auto const& w = ip_data.integration_weight;
- double T_int_pt = 0.0;
- double p_int_pt = 0.0;
- // Order matters: First T, then P!
- NumLib::shapeFunctionInterpolate(local_x, N, T_int_pt, p_int_pt);
+ double p_at_xi = 0.;
+ NumLib::shapeFunctionInterpolate(local_p, N, p_at_xi);
+ double T1_at_xi = 0.;
+ NumLib::shapeFunctionInterpolate(local_T1, N, T1_at_xi);
- // \todo the first argument has to be changed for non constant
- // porosity model
auto const porosity =
- process_data.porous_media_properties.getPorosity(t, pos).getValue(
- t, pos, 0.0, T_int_pt);
- auto const intrinsic_permeability =
- 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];
+ material_properties.porous_media_properties.getPorosity(t, pos)
+ .getValue(t, pos, 0.0, T1_at_xi);
vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)] =
- T_int_pt;
+ T1_at_xi;
vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::p)] =
- p_int_pt;
+ 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 =
- process_data.fluid_properties->getValue(
+ 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 thermal_dispersivity_longitudinal =
- process_data.thermal_dispersivity_longitudinal(t, pos)[0];
- auto const thermal_dispersivity_transversal =
- process_data.thermal_dispersivity_transversal(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 density = process_data.fluid_properties->getValue(
- MaterialLib::Fluid::FluidPropertyType::Density, vars);
+ 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;
+
+ local_M.noalias() += w * heat_capacity * 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 = process_data.fluid_properties->getValue(
+ auto const viscosity = material_properties.fluid_properties->getValue(
MaterialLib::Fluid::FluidPropertyType::Viscosity, vars);
- GlobalDimMatrixType K_over_mu = intrinsic_permeability / viscosity;
+ auto const intrinsic_permeability =
+ material_properties.porous_media_properties
+ .getIntrinsicPermeability(t, pos)
+ .getValue(t, pos, 0.0, T1_at_xi);
+
+ GlobalDimMatrixType K_over_mu = intrinsic_permeability / viscosity;
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);
+ material_properties.has_gravity
+ ? GlobalDimVectorType(-K_over_mu * (dNdx * local_p_Eigen_type -
+ fluid_density * b))
+ : GlobalDimVectorType(-K_over_mu * dNdx * local_p_Eigen_type);
double const velocity_magnitude = velocity.norm();
GlobalDimMatrixType const thermal_dispersivity =
@@ -163,36 +289,12 @@ void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
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
- Ktt.noalias() +=
- (dNdx.transpose() * hydrodynamic_thermodispersion * dNdx +
- N.transpose() * velocity.transpose() * dNdx *
- fluid_reference_density * specific_heat_capacity_fluid) *
- w;
- Kpp.noalias() += w * dNdx.transpose() * K_over_mu * dNdx;
- Mtt.noalias() += w * N.transpose() * heat_capacity * N;
- Mpp.noalias() += w * N.transpose() * specific_storage * N;
- if (process_data.has_gravity)
- Bp += w * density * dNdx.transpose() * K_over_mu * b;
- /* with Oberbeck-Boussing assumption density difference only exists
- * in buoyancy effects */
+ local_K.noalias() +=
+ w * (dNdx.transpose() * hydrodynamic_thermodispersion * dNdx +
+ N.transpose() * velocity.transpose() * dNdx * fluid_density *
+ specific_heat_capacity_fluid);
}
}
-template <typename ShapeFunction, typename IntegrationMethod,
- unsigned GlobalDim>
-void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
- assembleHeatTransportEquation(double const t,
- std::vector<double>& local_M_data,
- std::vector<double>& local_K_data,
- std::vector<double>& local_b_data,
- LocalCoupledSolutions const& coupled_term)
-{
-}
-
} // namespace HT
} // namespace ProcessLib
diff --git a/ProcessLib/HT/StaggeredHTFEM.h b/ProcessLib/HT/StaggeredHTFEM.h
index b6575cc5cb215c1ed533a633b9dfd7d367760e02..ca1217cee25e86bb2503aeb0ce96a51a3920d08c 100644
--- a/ProcessLib/HT/StaggeredHTFEM.h
+++ b/ProcessLib/HT/StaggeredHTFEM.h
@@ -14,7 +14,6 @@
#include <Eigen/Dense>
#include <vector>
-#include "HTProcessData.h"
#include "NumLib/DOF/DOFTableUtil.h"
#include "NumLib/Extrapolation/ExtrapolatableElement.h"
#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
@@ -56,10 +55,10 @@ public:
std::size_t const local_matrix_size,
bool is_axially_symmetric,
unsigned const integration_order,
- HTMaterialProperties const& process_data)
+ HTMaterialProperties const& material_properties)
: HTFEM<ShapeFunction, IntegrationMethod, GlobalDim>(
element, local_matrix_size, is_axially_symmetric,
- integration_order, process_data, 1)
+ integration_order, material_properties, 1)
{
}