Homec4science
Diffusion SpecMiCP / ReactMiCP (master) Browse

specmicp/tests/specmicp/adim/adimensional_system_sorption.cpp
856cb1763aadmaster

Last Change
Raw File

adimensional_system_sorption.cpp
View Options

#include "catch.hpp"
#include "specmicp_common/log.hpp"
#include "specmicp/adimensional/adimensional_system.hpp"
#include "specmicp_common/micpsolver/micpsolver.hpp"
#include "specmicp/adimensional/adimensional_system_solver.hpp"
#include "specmicp/problem_solver/formulation.hpp"
#include "specmicp/problem_solver/dissolver.hpp"
#include "specmicp/problem_solver/reactant_box.hpp"
#include "specmicp/adimensional/adimensional_system_solution_extractor.hpp"
#include "specmicp/io/adimensional_system_solution_saver.hpp"
#include "specmicp/io/adimensional_system_solution_reader.hpp"
#include "specmicp_database/database.hpp"
#include <iostream>
specmicp::RawDatabasePtr get_test_feoh_sorption_database()
{
specmicp::database::Database thedatabase(TEST_FEOH_PATH);
thedatabase.remove_half_cell_reactions();
return thedatabase.get_database();
}
static std::string ssites_appendix_database =
R"plop(
[
{
"label": ">SiOH",
}
]
)plop";
static std::string sorbed_appendix_database =
R"plop(
[
{
"label": ">SiO[-]",
"composition": "HO[-], -H2O, >SiOH",
"log_k": 12.70
},
{
"label": ">SiOCa[+]",
"composition": "Ca[2+], -H[+], >SiOH",
"log_k": 9.40
},
{
"label": ">SiOCaSO4[-]",
"composition": "Ca[2+], SO4[2-], -H[+], >SiOH",
"log_k": 6.00
},
{
"label": ">SiOCaCl",
"composition": "Ca[2+], Cl[-], -H[+], >SiOH",
"log_k": 8.90
},
{
"label": ">SiONa",
"composition": "Na[+], -H[+], >SiOH",
"log_k": 13.64
},
{
"label": ">SiOK",
"composition": "K[+], -H[+], >SiOH",
"log_k": 13.64
},
{
"label": ">SiOHCl[-]",
"composition": "Cl[-], >SiOH",
"log_k": 0.35
}
]
)plop";
specmicp::RawDatabasePtr get_test_cl_sorption_database()
{
specmicp::database::Database thedatabase(TEST_CEMDATA_PATH);
thedatabase.remove_half_cell_reactions();
std::map<std::string, std::string> swapping ({
{"H[+]","HO[-]"},
{"Si(OH)4", "SiO(OH)3[-]"},
{"Al[3+]","Al(OH)4[-]"},
});
thedatabase.swap_components(swapping);
thedatabase.add_surface_sites(ssites_appendix_database);
thedatabase.add_sorbed_species(sorbed_appendix_database);
thedatabase.minerals_keep_only({
"Portlandite",
"CSH,jennite",
"CSH,tobermorite",
"SiO2(am)",
"Calcite",
"Al(OH)3(mic)",
"C3AH6",
"C4AH19",
"C4AH13",
"C2AH7_5",
"CAH10",
"Monosulfoaluminate",
"Tricarboaluminate",
"Monocarboaluminate",
"Hemicarboaluminate",
"Gypsum",
"Ettringite",
"Friedels_salt"
});
return thedatabase.get_database();
}
using namespace specmicp;
TEST_CASE("Surface Sorption", "[sorption, specmicp, sorbed species, EDL]")
{
std::cerr.flush();
specmicp::logger::ErrFile::stream() = &std::cerr;
specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
SECTION("FeOH") {
RawDatabasePtr raw_data = get_test_feoh_sorption_database();
// index_t id_h2o = raw_data->get_id_component("H2O");
// index_t id_h = raw_data->get_id_component("H[+]");
// index_t id_hg = raw_data->get_id_component("Hg[2+]");
// index_t id_pb = raw_data->get_id_component("Pb[2+]");
// index_t id_so4 = raw_data->get_id_component("SO4[2-]");
// index_t id_na = raw_data->get_id_component("Na[+]");
// index_t id_cl = raw_data->get_id_component("Cl[-]");
units::UnitsSet tunits = units::SI_units;
tunits.length = units::LengthUnit::decimeter;
specmicp::ReactantBox reactant_box(raw_data, tunits);
reactant_box.set_solution(0.5, "kg/L");
//reactant_box.add_aqueous_species("NaOH", 1e-2, "mol/kg");
//reactant_box.add_aqueous_species("HCl", 1e-4, "mol/kg");
reactant_box.add_aqueous_species("NaCl", 0.01, "mol/kg");
reactant_box.add_aqueous_species("PbSO4", 100e-6, "mol/kg");
reactant_box.add_aqueous_species("HgSO4", 100e-6, "mol/kg");
reactant_box.add_component("H[+]");
Vector total_concentration = reactant_box.get_total_concentration(true);
std::cout << total_concentration << std::endl;
{
specmicp::Vector x;
specmicp::AdimensionalSystemConstraints constraints(total_concentration);
constraints.add_fixed_activity_component(raw_data->get_id_component("H[+]"), -4);
constraints.set_charge_keeper(raw_data->get_id_component("Cl[-]"));
Vector tot_surface (raw_data->nb_ssites());
tot_surface(raw_data->get_id_ssite(">Fe^wOH")) = 0.5*0.2*1.0/88.851;
tot_surface(raw_data->get_id_ssite(">Fe^sOH")) = 0.5*0.005*1.0/88.851;
//constraints.set_equilibrium_surface_model(tot_surface);
constraints.set_EDL_surface_model(tot_surface, 300e3);
//constraints.set_saturated_system();
specmicp::AdimensionalSystemSolver solver(raw_data, constraints);
solver.get_options().units_set = tunits;
solver.get_options().solver_options.max_iter = 200;
solver.get_options().solver_options.set_tolerance(1e-12);
solver.initialize_variables(x, 1.0, -6.0);
solver.get_options().units_set = tunits;
auto perf = solver.solve(x);
REQUIRE(perf.return_code >= micpsolver::MiCPSolverReturnCode::Success);
specmicp::AdimensionalSystemSolution sol = solver.get_raw_solution(x);
specmicp::AdimensionalSystemSolutionExtractor extr(sol, raw_data, tunits);
std::cout << " pH : " << extr.pH() << "\n ================== \n" << std::endl;
std::cout << "surface potential : " << extr.surface_potential()*1e3 << " mV " << std::endl;
std::cout << "surface charge : " << extr.surface_charge()*constants::faraday_constant/600*100 << " uC / cm^2 " << std::endl;
{
std::cout << " == weak == \n";
auto id_ws = raw_data->get_id_ssite(">Fe^wOH");
auto sum = 2*tot_surface(id_ws);//*extr.density_water();
//for (auto s: raw_data->range_sorbed()) {
// sum += raw_data->nu_ssites(s, id_ws)*sol.sorbed_molalities(s);
//}
std::cout << raw_data->get_label_ssite(id_ws) << ": " << extr.free_surface_concentration(id_ws) / sum * 100 << "\n";
for (auto s: raw_data->range_sorbed()) {
if (raw_data->nu_ssites(s, id_ws) != 0.0)
std::cout << raw_data->get_label_sorbed(s) << ": " << extr.molality_sorbed_species(s) / sum * 100 << "\n";
}
std::cout << std::endl;
}
{
std::cout << " == strong == \n";
auto id_ss = raw_data->get_id_ssite(">Fe^sOH");
auto sum = 2*tot_surface(id_ss);//*extr.density_water();
//for (auto s: raw_data->range_sorbed()) {
// sum += raw_data->nu_ssites(s, id_ws)*sol.sorbed_molalities(s);
//}
std::cout << raw_data->get_label_ssite(id_ss) << ": " << extr.free_surface_concentration(id_ss) / sum * 100 << "\n";
for (auto s: raw_data->range_sorbed()) {
if (raw_data->nu_ssites(s, id_ss) != 0.0)
std::cout << raw_data->get_label_sorbed(s) << ": " << extr.molality_sorbed_species(s) / sum * 100 << "\n";
}
std::cout << std::endl;
}
}
{
specmicp::Vector x;
specmicp::AdimensionalSystemConstraints constraints(total_concentration);
constraints.add_fixed_activity_component(raw_data->get_id_component("H[+]"), -8);
constraints.set_charge_keeper(raw_data->get_id_component("Cl[-]"));
Vector tot_surface (raw_data->nb_ssites());
tot_surface(raw_data->get_id_ssite(">Fe^wOH")) = 0.5*0.2*1.0/88.851;
tot_surface(raw_data->get_id_ssite(">Fe^sOH")) = 0.5*0.005*1.0/88.851;
//constraints.set_equilibrium_surface_model(tot_surface);
constraints.set_EDL_surface_model(tot_surface, 300e3);
//constraints.set_saturated_system();
specmicp::AdimensionalSystemSolver solver(raw_data, constraints);
solver.get_options().units_set = tunits;
solver.get_options().solver_options.max_iter = 200;
solver.get_options().solver_options.set_tolerance(1e-12);
solver.initialize_variables(x, 1.0, -6.0);
solver.get_options().units_set = tunits;
auto perf = solver.solve(x);
REQUIRE(perf.return_code >= micpsolver::MiCPSolverReturnCode::Success);
specmicp::AdimensionalSystemSolution sol = solver.get_raw_solution(x);
specmicp::AdimensionalSystemSolutionExtractor extr(sol, raw_data, tunits);
std::cout << " pH : " << extr.pH() << "\n ================== \n" << std::endl;
std::cout << "surface potential : " << extr.surface_potential()*1e3 << " mV " << std::endl;
std::cout << "surface charge : " << extr.surface_charge()*constants::faraday_constant/600*100 << " uC / cm^2 " << std::endl;
{
std::cout << " == weak == \n";
auto id_ws = raw_data->get_id_ssite(">Fe^wOH");
auto sum = 2*tot_surface(id_ws);///extr.density_water();
//for (auto s: raw_data->range_sorbed()) {
// sum += raw_data->nu_ssites(s, id_ws)*sol.sorbed_molalities(s);
//}
std::cout << raw_data->get_label_ssite(id_ws) << ": " << extr.free_surface_concentration(id_ws) / sum * 100 << "\n";
for (auto s: raw_data->range_sorbed()) {
if (raw_data->nu_ssites(s, id_ws) != 0.0)
std::cout << raw_data->get_label_sorbed(s) << ": " << extr.molality_sorbed_species(s) / sum * 100 << "\n";
}
std::cout << std::endl;
}
{
std::cout << " == strong == \n";
auto id_ss = raw_data->get_id_ssite(">Fe^sOH");
auto sum = 2*tot_surface(id_ss);///extr.density_water();
auto check_sum = 0.0;
//for (auto s: raw_data->range_sorbed()) {
// sum += raw_data->nu_ssites(s, id_ws)*sol.sorbed_molalities(s);
//}
std::cout << raw_data->get_label_ssite(id_ss) << ": " << extr.free_surface_concentration(id_ss) / sum * 100 << "\n";
check_sum += extr.free_surface_concentration(id_ss);
for (auto s: raw_data->range_sorbed()) {
if (raw_data->nu_ssites(s, id_ss) != 0.0) {
std::cout << raw_data->get_label_sorbed(s) << ": " << extr.molality_sorbed_species(s) / sum * 100 << "\n";
check_sum += extr.molality_sorbed_species(s);
}
}
std::cout << "Sum : " << check_sum/sum*100 << "\n";
std::cout << std::endl;
}
}
}
SECTION("NaCl") {
std::cout << "\n ------------------\n";
std::cout << " Cl sorption on CSH\n";
std::cout << " ------------------\n";
specmicp::RawDatabasePtr raw_data = get_test_cl_sorption_database();
specmicp::units::UnitsSet units_set;
units_set.length = specmicp::units::LengthUnit::decimeter;
//units_set.quantity = specmicp::units::QuantityUnit::millimoles;
specmicp::ReactantBox reactant_box(raw_data, units_set);
reactant_box.set_solution(0.25, "kg/L");
reactant_box.add_solid_phase("Portlandite", 4.46, "mol/L");
reactant_box.add_solid_phase("CSH,jennite", 4.59, "mol/L");
reactant_box.add_solid_phase("Monosulfoaluminate", 0.3, "mol/L");
reactant_box.add_solid_phase("Al(OH)3(mic)", 0.6, "mol/L");
reactant_box.add_solid_phase("Calcite", 0.2, "mol/L");
//reactant_box.add_aqueous_species("CO2", 0.5, "mol/L");
reactant_box.add_aqueous_species("KOH", 473.0, "mmol/L");
reactant_box.add_aqueous_species("NaOH", 52.0, "mmol/L");
reactant_box.add_component("Cl[-]");
reactant_box.set_fixed_saturation(0.7);
auto constraints = reactant_box.get_constraints(true);
specmicp::AdimensionalSystemSolver solver(raw_data, constraints);
solver.get_options().units_set = units_set;
solver.get_options().solver_options.set_maximum_iterations(200);
specmicp::Vector x;
solver.initialize_variables(x, 0.7, -6);
specmicp::micpsolver::MiCPPerformance perf = solver.solve(x, true);
REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet);
auto init_sol = solver.get_raw_solution(x);
AdimensionalSystemSolutionExtractor extr(init_sol, raw_data, units_set);
Vector sc(1);
sc(0) = 0.01*extr.mole_concentration_mineral(raw_data->get_id_mineral("CSH,jennite"));
constraints.set_EDL_surface_model(sc, 100e3);
solver = specmicp::AdimensionalSystemSolver(raw_data, constraints, init_sol);
solver.get_options().solver_options.set_tolerance(1e-10);
solver.get_options().units_set = units_set;
perf = solver.solve(x, true);
REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet);
auto sol = solver.get_raw_solution(x);
AdimensionalSystemSolutionExtractor extr_sorbed(sol, raw_data, units_set);
//
std::cout << " No NaCl \n ------ \n";
std::cout << " pH : " << extr.pH() << std::endl;
std::cout << "surface potential : " << extr_sorbed.surface_potential()*1e3 << " mV " << std::endl;
std::cout << "surface charge : " << extr_sorbed.surface_charge()*constants::faraday_constant/600*100 << " uC / cm^2 " << std::endl;
{
auto sum = sc(0)/(extr_sorbed.density_water()*extr_sorbed.porosity()*extr_sorbed.saturation_water());
std::cout << raw_data->get_label_ssite(0) << ": "
<< extr_sorbed.free_surface_concentration(0)
<< " = "
<< extr_sorbed.free_surface_concentration(0)/sum*100 << "\n";
for (auto s: raw_data->range_sorbed()) {
if (raw_data->nu_ssites(s, 0) != 0.0) {
std::cout << raw_data->get_label_sorbed(s)
<< ": " << extr_sorbed.molality_sorbed_species(s)
<< " = "
<< extr_sorbed.molality_sorbed_species(s)/sum*100
<< "\n";
}
}
}
constraints.add_fixed_molality_component(raw_data->get_id_component("Cl[-]"), std::log10(0.5));
constraints.set_charge_keeper(raw_data->get_id_component("Na[+]"));
solver = specmicp::AdimensionalSystemSolver(raw_data, constraints, sol);
solver.get_options().solver_options.set_tolerance(1e-10);
solver.get_options().units_set = units_set;
perf = solver.solve(x, true);
REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet);
auto sol2 = solver.get_raw_solution(x);
AdimensionalSystemSolutionExtractor extr_sorbed2(sol2, raw_data, units_set);
std::cout << "\n [Cl]=0.5M \n ------ \n";
std::cout << " pH : " << extr_sorbed2.pH() << std::endl;
std::cout << "surface potential : " << extr_sorbed2.surface_potential()*1e3 << " mV " << std::endl;
std::cout << "surface charge : " << extr_sorbed2.surface_charge()*constants::faraday_constant/600*100 << " uC / cm^2 " << std::endl;
{
auto sum = sc(0)/(extr_sorbed2.density_water()*extr_sorbed2.porosity()*extr_sorbed2.saturation_water());
std::cout << raw_data->get_label_ssite(0) << ": "
<< extr_sorbed2.free_surface_concentration(0)
<< " = "
<< extr_sorbed2.free_surface_concentration(0)/sum*100 << "\n";
for (auto s: raw_data->range_sorbed()) {
if (raw_data->nu_ssites(s, 0) != 0.0) {
std::cout << raw_data->get_label_sorbed(s)
<< ": " << extr_sorbed2.molality_sorbed_species(s)
<< " = "
<< extr_sorbed2.molality_sorbed_species(s)/sum*100
<< "\n";
}
}
}
}
}
c4science · Help