From 3069bb47dbc16f21307b6412745a786c5eb35349 Mon Sep 17 00:00:00 2001
From: renchao_lu <renchao.lu@gmail.com>
Date: Mon, 12 Apr 2021 17:32:08 +0200
Subject: [PATCH] [CL] Update volume fraction post reaction.
---
ChemistryLib/ChemicalSolverInterface.h | 8 ++++
ChemistryLib/PhreeqcIO.cpp | 38 +++++++++++++++++++
ChemistryLib/PhreeqcIO.h | 6 +++
.../ComponentTransportFEM.h | 24 +++++++++++-
4 files changed, 75 insertions(+), 1 deletion(-)
diff --git a/ChemistryLib/ChemicalSolverInterface.h b/ChemistryLib/ChemicalSolverInterface.h
index 728f24b8872..0f5533bc996 100644
--- a/ChemistryLib/ChemicalSolverInterface.h
+++ b/ChemistryLib/ChemicalSolverInterface.h
@@ -59,6 +59,14 @@ public:
return {};
}
+ virtual void updateVolumeFractionPostReaction(
+ GlobalIndexType const& /*chemical_system_id*/,
+ MaterialPropertyLib::Medium const* /*medium*/,
+ ParameterLib::SpatialPosition const& /*pos*/, double const /*porosity*/,
+ double const /*t*/, double const /*dt*/)
+ {
+ }
+
virtual void computeSecondaryVariable(
std::size_t const /*ele_id*/,
std::vector<GlobalIndexType> const& /*chemical_system_indices*/)
diff --git a/ChemistryLib/PhreeqcIO.cpp b/ChemistryLib/PhreeqcIO.cpp
index 86cd6e3e735..12a47ed160d 100644
--- a/ChemistryLib/PhreeqcIO.cpp
+++ b/ChemistryLib/PhreeqcIO.cpp
@@ -671,6 +671,44 @@ std::vector<std::string> const PhreeqcIO::getComponentList() const
return component_names;
}
+void PhreeqcIO::updateVolumeFractionPostReaction(
+ GlobalIndexType const& chemical_system_id,
+ MaterialPropertyLib::Medium const* medium,
+ ParameterLib::SpatialPosition const& pos, double const porosity,
+ double const t, double const dt)
+{
+ auto const& solid_phase = medium->phase("Solid");
+ auto const& liquid_phase = medium->phase("AqueousLiquid");
+
+ MaterialPropertyLib::VariableArray vars;
+
+ auto const liquid_density =
+ liquid_phase.property(MaterialPropertyLib::PropertyType::density)
+ .template value<double>(vars, pos, t, dt);
+
+ for (auto& kinetic_reactant : _chemical_system->kinetic_reactants)
+ {
+ auto const& solid_constituent =
+ solid_phase.component(kinetic_reactant.name);
+
+ if (solid_constituent.hasProperty(
+ MaterialPropertyLib::PropertyType::molality))
+ {
+ continue;
+ }
+
+ auto const molar_volume =
+ solid_constituent
+ .property(MaterialPropertyLib::PropertyType::molar_volume)
+ .template value<double>(vars, pos, t, dt);
+
+ (*kinetic_reactant.volume_fraction)[chemical_system_id] +=
+ ((*kinetic_reactant.molality)[chemical_system_id] -
+ (*kinetic_reactant.molality_prev)[chemical_system_id]) *
+ liquid_density * porosity * molar_volume;
+ }
+}
+
void PhreeqcIO::computeSecondaryVariable(
std::size_t const ele_id,
std::vector<GlobalIndexType> const& chemical_system_indices)
diff --git a/ChemistryLib/PhreeqcIO.h b/ChemistryLib/PhreeqcIO.h
index e78b048a971..01bafd23313 100644
--- a/ChemistryLib/PhreeqcIO.h
+++ b/ChemistryLib/PhreeqcIO.h
@@ -72,6 +72,12 @@ public:
friend std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io);
+ void updateVolumeFractionPostReaction(
+ GlobalIndexType const& chemical_system_id,
+ MaterialPropertyLib::Medium const* medium,
+ ParameterLib::SpatialPosition const& pos, double const porosity,
+ double const t, double const dt) override;
+
void computeSecondaryVariable(
std::size_t const ele_id,
std::vector<GlobalIndexType> const& chemical_system_indices) override;
diff --git a/ProcessLib/ComponentTransport/ComponentTransportFEM.h b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
index 95d1b686b40..52b1d800a2b 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportFEM.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
@@ -1158,7 +1158,7 @@ public:
void computeSecondaryVariableConcrete(
double const t,
- double const /*dt*/,
+ double const dt,
Eigen::VectorXd const& local_x,
Eigen::VectorXd const& /*local_x_dot*/) override
{
@@ -1183,6 +1183,28 @@ public:
if (_process_data.chemical_solver_interface)
{
+ if (_process_data.chemically_induced_porosity_change)
+ {
+ auto const& medium = _process_data.media_map->getMedium(ele_id);
+
+ ParameterLib::SpatialPosition pos;
+ pos.setElementID(ele_id);
+
+ for (auto& ip_data : _ip_data)
+ {
+ ip_data.porosity = ip_data.porosity_prev;
+
+ _process_data.chemical_solver_interface
+ ->updateVolumeFractionPostReaction(
+ ip_data.chemical_system_id, medium, pos,
+ ip_data.porosity, t, dt);
+
+ _process_data.chemical_solver_interface
+ ->updatePorosityPostReaction(chemical_system_id, medium,
+ porosity);
+ }
+ }
+
std::vector<GlobalIndexType> chemical_system_indices;
chemical_system_indices.reserve(n_integration_points);
std::transform(
--
GitLab