diff --git a/ChemistryLib/PhreeqcIO.cpp b/ChemistryLib/PhreeqcIO.cpp index a93dda604377772d05b0a7f273d14452be3d30de..a4228a09c0319015e22b73a06f5ab55a0baf277e 100644 --- a/ChemistryLib/PhreeqcIO.cpp +++ b/ChemistryLib/PhreeqcIO.cpp @@ -703,6 +703,28 @@ void PhreeqcIO::updateVolumeFractionPostReaction( (*kinetic_reactant.molality_prev)[chemical_system_id]) * liquid_density * porosity * molar_volume; } + + for (auto& equilibrium_reactant : _chemical_system->equilibrium_reactants) + { + auto const& solid_constituent = + solid_phase.component(equilibrium_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); + + (*equilibrium_reactant.volume_fraction)[chemical_system_id] += + ((*equilibrium_reactant.molality)[chemical_system_id] - + (*equilibrium_reactant.molality_prev)[chemical_system_id]) * + liquid_density * porosity * molar_volume; + } } void PhreeqcIO::updatePorosityPostReaction( @@ -726,6 +748,22 @@ void PhreeqcIO::updatePorosityPostReaction( ((*kinetic_reactant.volume_fraction)[chemical_system_id] - (*kinetic_reactant.volume_fraction_prev)[chemical_system_id]); } + + for (auto& equilibrium_reactant : _chemical_system->equilibrium_reactants) + { + auto const& solid_constituent = + solid_phase.component(equilibrium_reactant.name); + + if (solid_constituent.hasProperty( + MaterialPropertyLib::PropertyType::molality)) + { + continue; + } + + porosity -= + ((*equilibrium_reactant.volume_fraction)[chemical_system_id] - + (*equilibrium_reactant.volume_fraction_prev)[chemical_system_id]); + } } void PhreeqcIO::computeSecondaryVariable(