diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp
index d3d46f9b44269b2af7d50493a89857026d2a860f..e1815115ba087dc573fc271ac26ef4d0c2ccdb50 100644
--- a/Applications/ApplicationsLib/ProjectData.cpp
+++ b/Applications/ApplicationsLib/ProjectData.cpp
@@ -35,6 +35,7 @@
#include "ProcessLib/UncoupledProcessesTimeLoop.h"
#include "ProcessLib/GroundwaterFlow/CreateGroundwaterFlowProcess.h"
+#include "ProcessLib/HC/CreateHCProcess.h"
#include "ProcessLib/HT/CreateHTProcess.h"
#include "ProcessLib/HeatConduction/CreateHeatConductionProcess.h"
#include "ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.h"
@@ -383,6 +384,13 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config,
_process_variables, _parameters, integration_order,
process_config);
}
+ else if (type == "HC")
+ {
+ process = ProcessLib::HC::createHCProcess(
+ *_mesh_vec[0], std::move(jacobian_assembler),
+ _process_variables, _parameters, integration_order,
+ process_config);
+ }
else if (type == "SMALL_DEFORMATION")
{
//! \ogs_file_param{prj__processes__process__SMALL_DEFORMATION__dimension}
diff --git a/ProcessLib/CMakeLists.txt b/ProcessLib/CMakeLists.txt
index befb3595358128e00fdeb2e5732c9ab60572756e..9e66da45efbff37523ee5ed9ee8eb168d9287b73 100644
--- a/ProcessLib/CMakeLists.txt
+++ b/ProcessLib/CMakeLists.txt
@@ -45,6 +45,9 @@ APPEND_SOURCE_FILES(SOURCES HeatConduction)
add_subdirectory(RichardsFlow)
APPEND_SOURCE_FILES(SOURCES RichardsFlow)
+add_subdirectory(HC)
+APPEND_SOURCE_FILES(SOURCES HC)
+
add_subdirectory(HT)
APPEND_SOURCE_FILES(SOURCES HT)
diff --git a/ProcessLib/HC/CreateHCProcess.cpp b/ProcessLib/HC/CreateHCProcess.cpp
index 383c897063b367b99e823daf258468544115dcc3..078fcfee26c17ef17f6b8435fde6be75f3e6a831 100644
--- a/ProcessLib/HC/CreateHCProcess.cpp
+++ b/ProcessLib/HC/CreateHCProcess.cpp
@@ -63,13 +63,6 @@ std::unique_ptr<Process> createHCProcess(
auto viscosity_model =
MaterialLib::Fluid::createViscosityModel(viscosity_conf);
- // Parameter for the density of the solid.
- auto& density_solid = findParameter<double>(
- config,
- //! \ogs_file_param_special{prj__processes__process__HC__density_solid}
- "density_solid", parameters, 1);
- DBUG("Use \'%s\' as density_solid parameter.", density_solid.name.c_str());
-
//! \ogs_file_param{prj__processes__process__HC__fluid__density}
auto const& fluid_density_conf = fluid_config.getConfigSubtree("density");
auto fluid_density =
@@ -83,56 +76,44 @@ std::unique_ptr<Process> createHCProcess(
DBUG("Use \'%s\' as fluid_reference_density parameter.",
fluid_reference_density.name.c_str());
- // Parameter for the specific heat capacity of the solid.
- auto& specific_heat_capacity_solid = findParameter<double>(
- config,
- //! \ogs_file_param_special{prj__processes__process__HC__specific_heat_capacity_solid}
- "specific_heat_capacity_solid", parameters, 1);
- DBUG("Use \'%s\' as specific_heat_capacity_solid parameter.",
- specific_heat_capacity_solid.name.c_str());
-
- // Parameter for the specific heat capacity of the fluid.
- auto& specific_heat_capacity_fluid = findParameter<double>(
+ // Parameter for the longitudinal solute dispersivity.
+ auto const& molecular_diffusion_coefficient = findParameter<double>(
config,
- //! \ogs_file_param_special{prj__processes__process__HC__specific_heat_capacity_fluid}
- "specific_heat_capacity_fluid", parameters, 1);
- DBUG("Use \'%s\' as specific_heat_capacity_fluid parameter.",
- specific_heat_capacity_fluid.name.c_str());
-
- // Parameter for the thermal conductivity of the solid (only one scalar per
- // element, i.e., the isotropic case is handled at the moment)
- auto& thermal_dispersivity_longitudinal = findParameter<double>(
+ //!
+ //\ogs_file_param_special{prj__processes__process__HC__molecular_diffusion_coefficient
+ "molecular_diffusion_coefficient", parameters, 1);
+ DBUG("Use \'%s\' as molecular diffusion coefficient parameter.",
+ molecular_diffusion_coefficient.name.c_str());
+
+ // Parameter for the longitudinal solute dispersivity.
+ auto const& solute_dispersivity_longitudinal = findParameter<double>(
config,
- //! \ogs_file_param_special{prj__processes__process__HC__thermal_dispersivity_longitudinal}
- "thermal_dispersivity_longitudinal", parameters, 1);
- DBUG("Use \'%s\' as thermal_dispersivity_longitudinal parameter.",
- thermal_dispersivity_longitudinal.name.c_str());
-
- // Parameter for the thermal conductivity of the solid (only one scalar per
- // element, i.e., the isotropic case is handled at the moment)
- auto& thermal_dispersivity_transversal = findParameter<double>(
- config,
- //! \ogs_file_param_special{prj__processes__process__HC__thermal_dispersivity_transversal}
- "thermal_dispersivity_transversal", parameters, 1);
- DBUG("Use \'%s\' as thermal_dispersivity_transversal parameter.",
- thermal_dispersivity_transversal.name.c_str());
-
- // Parameter for the thermal conductivity of the solid (only one scalar per
- // element, i.e., the isotropic case is handled at the moment)
- auto& thermal_conductivity_solid = findParameter<double>(
- config,
- //! \ogs_file_param_special{prj__processes__process__HC__thermal_conductivity_solid}
- "thermal_conductivity_solid", parameters, 1);
- DBUG("Use \'%s\' as thermal_conductivity_solid parameter.",
- thermal_conductivity_solid.name.c_str());
-
- // Parameter for the thermal conductivity of the fluid.
- auto& thermal_conductivity_fluid = findParameter<double>(
+ //!
+ //\ogs_file_param_special{prj__processes__process__HC__solute_dispersivity_longitudinal
+ "solute_dispersivity_longitudinal", parameters, 1);
+ DBUG("Use \'%s\' as longitudinal solute dispersivity parameter.",
+ solute_dispersivity_longitudinal.name.c_str());
+
+ // Parameter for the transverse solute dispersivity.
+ auto const& solute_dispersivity_transverse = findParameter<double>(
config,
- //! \ogs_file_param_special{prj__processes__process__HC__thermal_conductivity_fluid}
- "thermal_conductivity_fluid", parameters, 1);
- DBUG("Use \'%s\' as thermal_conductivity_fluid parameter.",
- thermal_conductivity_fluid.name.c_str());
+ //!
+ //\ogs_file_param_special{prj__processes__process__HC__solute_dispersivity_transverse
+ "solute_dispersivity_transverse", parameters, 1);
+ DBUG("Use \'%s\' as transverse solute dispersivity parameter.",
+ solute_dispersivity_transverse.name.c_str());
+
+ // Parameter for the retardation factor.
+ auto const& retardation_factor =
+ findParameter<double>(config,
+ //! \ogs_file_param_special{prj__processes__process__HC__retardation_factor}
+ "retardation_factor", parameters, 1);
+
+ // Parameter for the decay rate.
+ auto const& decay_rate =
+ findParameter<double>(config,
+ //! \ogs_file_param_special{prj__processes__process__HC__decay_rate}
+ "decay_rate", parameters, 1);
// Specific body force parameter.
Eigen::Vector3d specific_body_force;
@@ -147,22 +128,20 @@ std::unique_ptr<Process> createHCProcess(
HCProcessData process_data{
std::move(porous_media_properties),
std::move(viscosity_model),
- density_solid,
fluid_reference_density,
std::move(fluid_density),
- thermal_dispersivity_longitudinal,
- thermal_dispersivity_transversal,
- specific_heat_capacity_solid,
- specific_heat_capacity_fluid,
- thermal_conductivity_solid,
- thermal_conductivity_fluid,
+ molecular_diffusion_coefficient,
+ solute_dispersivity_longitudinal,
+ solute_dispersivity_transverse,
+ retardation_factor,
+ decay_rate,
specific_body_force,
has_gravity};
SecondaryVariableCollection secondary_variables;
NumLib::NamedFunctionCaller named_function_caller(
- {"HC_temperature_pressure"});
+ {"HC_concentration_pressure"});
ProcessLib::parseSecondaryVariables(config, secondary_variables,
named_function_caller);
diff --git a/ProcessLib/HC/HCFEM.h b/ProcessLib/HC/HCFEM.h
index ede8a2e9c5a837816b00b7bf306fd868ef6b2dfe..6464581e7b970f369f8873392f6cd684aa01edbb 100644
--- a/ProcessLib/HC/HCFEM.h
+++ b/ProcessLib/HC/HCFEM.h
@@ -27,7 +27,7 @@ namespace ProcessLib
{
namespace HC
{
-template < typename NodalRowVectorType, typename GlobalDimNodalMatrixType>
+template <typename NodalRowVectorType, typename GlobalDimNodalMatrixType>
struct IntegrationPointData final
{
IntegrationPointData(NodalRowVectorType const& N_,
@@ -148,70 +148,59 @@ public:
MaterialLib::Fluid::FluidProperty::ArrayType vars;
- for (std::size_t ip(0); ip < n_integration_points; ip++)
+ auto KCC = local_K.template block<num_nodes, num_nodes>(0, 0);
+ auto MCC = local_M.template block<num_nodes, num_nodes>(0, 0);
+ auto Kpp =
+ local_K.template block<num_nodes, num_nodes>(num_nodes, num_nodes);
+ auto Mpp =
+ local_M.template block<num_nodes, num_nodes>(num_nodes, num_nodes);
+ auto Bp = local_b.template block<num_nodes, 1>(num_nodes, 0);
+
+ for (std::size_t ip(0); ip < n_integration_points; ++ip)
{
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 = _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 C_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);
+ // Order matters: First C, then p!
+ NumLib::shapeFunctionInterpolate(local_x, N, C_int_pt, p_int_pt);
// \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(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];
- auto const specific_heat_capacity_fluid =
- _process_data.specific_heat_capacity_fluid(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 Ktt = local_K.template block<num_nodes, num_nodes>(0, 0);
- auto Mtt = local_M.template block<num_nodes, num_nodes>(0, 0);
- auto Kpp = local_K.template block<num_nodes, num_nodes>(num_nodes,
- num_nodes);
- auto Mpp = local_M.template block<num_nodes, num_nodes>(num_nodes,
- num_nodes);
- auto Bp = local_b.template block<num_nodes, 1>(num_nodes, 0);
+ .getValue(0.0, C_int_pt);
+
+ auto const retardation_factor =
+ _process_data.retardation_factor(t, pos)[0];
+
+
+ auto const& solute_dispersivity_transverse =
+ _process_data.solute_dispersivity_transverse(t, pos)[0];
+ auto const& solute_dispersivity_longitudinal =
+ _process_data.solute_dispersivity_longitudinal(t, pos)[0];
// Use the fluid density model to compute the density
vars[static_cast<int>(
- MaterialLib::Fluid::PropertyVariableType::T)] = T_int_pt;
+ MaterialLib::Fluid::PropertyVariableType::C)] = C_int_pt;
vars[static_cast<int>(
MaterialLib::Fluid::PropertyVariableType::p)] = p_int_pt;
- auto const density_water_T =
+ auto const density_water =
_process_data.fluid_density->getValue(vars);
+ auto const& decay_rate = _process_data.decay_rate(t, pos)[0];
+ auto const& molecular_diffusion_coefficient =
+ _process_data.molecular_diffusion_coefficient(t, pos)[0];
auto const& K =
_process_data.porous_media_properties.getIntrinsicPermeability(
@@ -223,44 +212,42 @@ public:
GlobalDimMatrixType const K_over_mu = K / mu;
GlobalDimVectorType const velocity =
- -perm_over_visc * (dNdx * p_nodal_values - density_water_T * b);
+ -perm_over_visc * (dNdx * p_nodal_values - density_water * b);
double const velocity_magnitude = velocity.norm();
GlobalDimMatrixType const& I(
GlobalDimMatrixType::Identity(GlobalDim, GlobalDim));
- GlobalDimMatrixType 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());
-
- auto 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;
+ GlobalDimMatrixType const hydrodynamic_dispersion =
+ velocity_magnitude != 0.0
+ ? GlobalDimMatrixType(
+ (porosity * molecular_diffusion_coefficient +
+ solute_dispersivity_transverse *
+ velocity_magnitude) *
+ I +
+ (solute_dispersivity_longitudinal -
+ solute_dispersivity_transverse) /
+ velocity_magnitude * velocity *
+ velocity.transpose())
+ : GlobalDimMatrixType(
+ (porosity * molecular_diffusion_coefficient +
+ solute_dispersivity_transverse *
+ velocity_magnitude) *
+ I);
// matrix assembly
- Ktt.noalias() +=
- (dNdx.transpose() * hydrodynamic_thermodispersion * dNdx +
- N.transpose() * velocity.transpose() * dNdx *
- fluid_reference_density * specific_heat_capacity_fluid) *
+ KCC.noalias() +=
+ (dNdx.transpose() * hydrodynamic_dispersion * dNdx +
+ N.transpose() * velocity.transpose() * dNdx +
+ N.transpose() * decay_rate * porosity * retardation_factor *
+ N) *
w;
+ MCC.noalias() +=
+ w * N.transpose() * porosity * retardation_factor * N;
Kpp.noalias() += w * dNdx.transpose() * perm_over_visc * dNdx;
- Mtt.noalias() += w * N.transpose() * heat_capacity * N;
Mpp.noalias() += w * N.transpose() * specific_storage * N;
- Bp += w * density_water_T * dNdx.transpose() * perm_over_visc * b;
+ Bp += w * density_water * dNdx.transpose() * perm_over_visc * b;
/* with Oberbeck-Boussing assumption density difference only exists
* in buoyancy effects */
-
- // velocity computed for output.
- for (unsigned d = 0; d < GlobalDim; ++d)
- {
- _darcy_velocities[d][ip] = velocity[d];
- }
}
}
diff --git a/ProcessLib/HC/HCProcess.cpp b/ProcessLib/HC/HCProcess.cpp
index 9323747685e334ae10af00b2712b168d2bc8a8cd..cb84dc12efb98e9b73314e6a59775b9f5947d1bf 100644
--- a/ProcessLib/HC/HCProcess.cpp
+++ b/ProcessLib/HC/HCProcess.cpp
@@ -67,13 +67,13 @@ void HCProcess::initializeConcreteProcess(
}
}
-void HCProcess::assembleConcreteProcess(const double t,
- GlobalVector const& x,
- GlobalMatrix& M,
- GlobalMatrix& K,
- GlobalVector& b,
- StaggeredCouplingTerm const&
- coupling_term)
+void HCProcess::assembleConcreteProcess(
+ const double t,
+ GlobalVector const& x,
+ GlobalMatrix& M,
+ GlobalMatrix& K,
+ GlobalVector& b,
+ StaggeredCouplingTerm const& coupling_term)
{
DBUG("Assemble HCProcess.");
diff --git a/ProcessLib/HC/HCProcess.h b/ProcessLib/HC/HCProcess.h
index 2a6dbd3e01c8e35d4732d13f8d34c974b73132f0..4e72ad050a5fd968e384a7658c997293ac3d2ddc 100644
--- a/ProcessLib/HC/HCProcess.h
+++ b/ProcessLib/HC/HCProcess.h
@@ -63,11 +63,9 @@ private:
MeshLib::Mesh const& mesh,
unsigned const integration_order) override;
- void assembleConcreteProcess(const double t, GlobalVector const& x,
- GlobalMatrix& M, GlobalMatrix& K,
- GlobalVector& b,
- StaggeredCouplingTerm const& coupling_term
- ) override;
+ void assembleConcreteProcess(
+ const double t, GlobalVector const& x, GlobalMatrix& M, GlobalMatrix& K,
+ GlobalVector& b, StaggeredCouplingTerm const& coupling_term) override;
void assembleWithJacobianConcreteProcess(
const double t, GlobalVector const& x, GlobalVector const& xdot,
diff --git a/ProcessLib/HC/HCProcessData.h b/ProcessLib/HC/HCProcessData.h
index 09efb4039ce898e1b7a79a945a04bf421da44872..d92119210ffdc97b2eb6dc550b8e4c0c714159ee 100644
--- a/ProcessLib/HC/HCProcessData.h
+++ b/ProcessLib/HC/HCProcessData.h
@@ -34,28 +34,24 @@ struct HCProcessData
HCProcessData(
PorousMediaProperties&& porous_media_properties_,
std::unique_ptr<MaterialLib::Fluid::FluidProperty>&& viscosity_model_,
- ProcessLib::Parameter<double> const& density_solid_,
ProcessLib::Parameter<double> const& fluid_reference_density_,
std::unique_ptr<MaterialLib::Fluid::FluidProperty>&& fluid_density_,
- ProcessLib::Parameter<double> const& thermal_dispersivity_longitudinal_,
- ProcessLib::Parameter<double> const& thermal_dispersivity_transversal_,
- ProcessLib::Parameter<double> const& specific_heat_capacity_solid_,
- ProcessLib::Parameter<double> const& specific_heat_capacity_fluid_,
- ProcessLib::Parameter<double> const& thermal_conductivity_solid_,
- ProcessLib::Parameter<double> const& thermal_conductivity_fluid_,
+ ProcessLib::Parameter<double> const& molecular_diffusion_coefficient_,
+ ProcessLib::Parameter<double> const& solute_dispersivity_longitudinal_,
+ ProcessLib::Parameter<double> const& solute_dispersivity_transverse_,
+ ProcessLib::Parameter<double> const& retardation_factor_,
+ ProcessLib::Parameter<double> const& decay_rate_,
Eigen::Vector3d const& specific_body_force_,
bool const has_gravity_)
: porous_media_properties(std::move(porous_media_properties_)),
viscosity_model(std::move(viscosity_model_)),
- density_solid(density_solid_),
fluid_reference_density(fluid_reference_density_),
fluid_density(std::move(fluid_density_)),
- specific_heat_capacity_solid(specific_heat_capacity_solid_),
- specific_heat_capacity_fluid(specific_heat_capacity_fluid_),
- thermal_dispersivity_longitudinal(thermal_dispersivity_longitudinal_),
- thermal_dispersivity_transversal(thermal_dispersivity_transversal_),
- thermal_conductivity_solid(thermal_conductivity_solid_),
- thermal_conductivity_fluid(thermal_conductivity_fluid_),
+ molecular_diffusion_coefficient(molecular_diffusion_coefficient_),
+ solute_dispersivity_longitudinal(solute_dispersivity_longitudinal_),
+ solute_dispersivity_transverse(solute_dispersivity_transverse_),
+ retardation_factor(retardation_factor_),
+ decay_rate(decay_rate_),
specific_body_force(specific_body_force_),
has_gravity(has_gravity_)
{
@@ -64,17 +60,15 @@ struct HCProcessData
HCProcessData(HCProcessData&& other)
: porous_media_properties(std::move(other.porous_media_properties)),
viscosity_model(other.viscosity_model.release()),
- density_solid(other.density_solid),
fluid_reference_density(other.fluid_reference_density),
fluid_density(other.fluid_density.release()),
- specific_heat_capacity_solid(other.specific_heat_capacity_solid),
- specific_heat_capacity_fluid(other.specific_heat_capacity_fluid),
- 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),
+ molecular_diffusion_coefficient(
+ other.molecular_diffusion_coefficient),
+ solute_dispersivity_longitudinal(
+ other.solute_dispersivity_longitudinal),
+ solute_dispersivity_transverse(other.solute_dispersivity_transverse),
+ retardation_factor(other.retardation_factor),
+ decay_rate(other.decay_rate),
specific_body_force(other.specific_body_force),
has_gravity(other.has_gravity)
{
@@ -91,15 +85,13 @@ struct HCProcessData
PorousMediaProperties porous_media_properties;
std::unique_ptr<MaterialLib::Fluid::FluidProperty> viscosity_model;
- Parameter<double> const& density_solid;
Parameter<double> const& fluid_reference_density;
std::unique_ptr<MaterialLib::Fluid::FluidProperty> fluid_density;
- Parameter<double> const& specific_heat_capacity_solid;
- Parameter<double> const& specific_heat_capacity_fluid;
- Parameter<double> const& thermal_dispersivity_longitudinal;
- Parameter<double> const& thermal_dispersivity_transversal;
- Parameter<double> const& thermal_conductivity_solid;
- Parameter<double> const& thermal_conductivity_fluid;
+ Parameter<double> const& molecular_diffusion_coefficient;
+ Parameter<double> const& solute_dispersivity_longitudinal;
+ Parameter<double> const& solute_dispersivity_transverse;
+ Parameter<double> const& retardation_factor;
+ Parameter<double> const& decay_rate;
Eigen::Vector3d const specific_body_force;
bool const has_gravity;
};