Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F90904996
acc_carbo_variables.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Tue, Nov 5, 20:18
Size
7 KB
Mime Type
text/x-c++
Expires
Thu, Nov 7, 20:18 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
22156362
Attached To
rSPECMICP SpecMiCP / ReactMiCP
acc_carbo_variables.cpp
View Options
#include "reactmicp_unsaturated.hpp"
#include "physics/unsaturated_laws.hpp"
#include "utils/plugins/plugin_interface.h"
#include "utils/plugins/plugin_base.hpp"
#include "utils/compat.hpp"
#include <Eigen/QR>
using namespace specmicp;
using namespace specmicp::reactmicp::solver;
using namespace specmicp::reactmicp::systems::unsaturated;
AdimensionalSystemSolution get_initial_condition(
database::RawDatabasePtr the_database,
const units::UnitsSet& units_set
);
void compo_from_oxyde(Vector& compo_oxyde, Vector& compo_species);
scalar_t vapor_pressure(index_t node, scalar_t saturation);
class SPECMICP_DLL_PUBLIC AccCarboVariablesFactory: public InitializeVariables
{
public:
static AccCarboVariablesFactory* get() {
return new AccCarboVariablesFactory();
}
virtual std::shared_ptr<UnsaturatedVariables> initialize_variables(
std::shared_ptr<database::DataContainer> the_database,
std::shared_ptr<mesh::Mesh1D> the_mesh,
const units::UnitsSet& the_units,
io::YAMLConfigHandle& configuration
) override
{
AdimensionalSystemSolution init = get_initial_condition(
the_database,
the_units);
AdimensionalSystemSolutionExtractor extr(init, the_database, the_units);
index_t id_co2 = the_database->get_id_component("CO2");
VariablesInterface variable_interface(the_mesh, the_database, {0, id_co2}, the_units);
variable_interface.set_binary_gas_diffusivity(0, 0.240);
variable_interface.set_binary_gas_diffusivity(id_co2, 0.160);
variable_interface.set_vapor_pressure_model(&vapor_pressure);
for (index_t node=1; node<the_mesh->nb_nodes(); ++node)
{
variable_interface.initialize_variables(node, extr);
}
variable_interface.set_porosity(0, 1.0);
variable_interface.set_liquid_saturation(0, 0.0);
variable_interface.set_partial_pressure(id_co2, 0, 1000);
variable_interface.set_partial_pressure(0, 0, 1000);
std::shared_ptr<UnsaturatedVariables> vars = variable_interface.get_variables();
vars->set_aqueous_scaling(0, 1.0e3/18.3e-3*1e-6);
for (auto aqc: the_database->range_aqueous_component())
{
vars->set_aqueous_scaling(aqc, 100e-6);
}
return vars;
}
};
class SPECMICP_DLL_PUBLIC AccCarboVariablesPlugin: public specmicp::plugins::PluginBase
{
public:
AccCarboVariablesPlugin():
PluginBase()
{
set_api_version(0, 1, 0);
}
bool initialize(const specmicp::plugins::PluginManagerServices& services) override
{
auto retcode = services.register_object(
"variables_factory", "acc_carbo_variables",
&AccCarboVariablesFactory::get
);
if (not retcode) {
return false;
}
return true;
}
};
SPECMICP_PLUGIN(AccCarboVariablesPlugin)
// implementation
scalar_t capillary_pressure(index_t node, scalar_t saturation)
{
if (saturation == 0.0) return 0.0;
return laws::capillary_pressure_BaroghelBouny(saturation, 37.5479e6, 2.1684);
}
scalar_t vapor_pressure(index_t node, scalar_t saturation)
{
scalar_t pv;
if (saturation == 0.0) pv = 0.0;
else if (saturation >= 1.0) pv = 3170;
else pv = 3170*laws::invert_kelvin_equation(
capillary_pressure(node, saturation));
return pv;
}
AdimensionalSystemSolution get_initial_condition(
database::RawDatabasePtr the_database,
const units::UnitsSet& units_set
)
{
Vector oxyde_compo(4);
oxyde_compo << 67.98, 22.66, 5.19, 2.96;
Vector species_compo;
compo_from_oxyde(oxyde_compo, species_compo);
scalar_t mult = 1e-6;
//species_compo *= 0.947e-2; // poro 0.47
//species_compo *= 1.08e-2; // poro 0.4
species_compo *= 1.25e-2; // poro 0.3
//species_compo *= 1.1e-2;
Formulation formulation;
formulation.mass_solution = 5e-4;
formulation.amount_minerals = {
{"C3S", species_compo(0)},
{"C2S", species_compo(1)},
{"C3A", species_compo(2)},
{"Gypsum", species_compo(3)},
{"Calcite", species_compo(0)*1e-3}
};
//formulation.extra_components_to_keep = {"HCO3[-]"};
formulation.concentration_aqueous =
{
{"NaOH", mult*0.0298},
{"KOH", mult*0.0801},
{"Cl[-]", mult*0.0001}
};
Dissolver dissolver(the_database);
AdimensionalSystemConstraints constraints;
constraints.set_total_concentrations(
dissolver.dissolve(formulation, true)
);
//constraints.set_inert_volume_fraction(0.1);
AdimensionalSystemSolver the_solver(the_database, constraints);
the_solver.get_options().units_set = units_set;
the_solver.get_options().system_options.non_ideality_tolerance = 1e-12;
micpsolver::MiCPSolverOptions* solver_options = &the_solver.get_options().solver_options;
//solver_options.maxstep = 10.0;
solver_options->maxiter_maxstep = 200;
solver_options->max_iter = 200;
solver_options->set_tolerance(1e-9);
Vector x;
the_solver.initialise_variables(x, 0.8, -4);
micpsolver::MiCPPerformance perf = the_solver.solve(x);
if (perf.return_code < micpsolver::MiCPSolverReturnCode::Success)
{
throw std::runtime_error("Error : failed to solve initial conditions");
}
auto prev_sol = the_solver.get_raw_solution(x);
prev_sol.main_variables(0) *= 0.8;
constraints.set_water_partial_pressure_model(
std::bind(&vapor_pressure,
1,
std::placeholders::_1)
);
the_solver = AdimensionalSystemSolver (the_database, constraints, prev_sol);
the_solver.get_options().units_set = units_set;
the_solver.get_options().system_options.non_ideality_tolerance = 1e-12;
solver_options = &the_solver.get_options().solver_options;
//solver_options.maxstep = 10.0;
solver_options->maxiter_maxstep = 200;
solver_options->set_tolerance(1e-9);
perf = the_solver.solve(x);
if (perf.return_code < micpsolver::MiCPSolverReturnCode::Success)
{
throw std::runtime_error("Error : failed to solve initial conditions");
}
return the_solver.get_raw_solution(x);
}
void compo_from_oxyde(Vector& compo_oxyde, Vector& compo_species)
{
constexpr double M_CaO = 56.08;
constexpr double M_SiO2 = 60.09;
constexpr double M_Al2O3 = 101.96;
constexpr double M_SO3 = 80.06;
Eigen::MatrixXd compo_in_oxyde(4, 4);
Vector molar_mass(4);
molar_mass << 1.0/M_CaO, 1.0/M_SiO2, 1.0/M_Al2O3, 1.0/M_SO3;
compo_oxyde = compo_oxyde.cwiseProduct(molar_mass);
//C3S C2S C3A Gypsum
compo_in_oxyde << 3, 2, 3, 1, // CaO
1, 1, 0, 0, // SiO2
0, 0, 1, 0, // Al2O3
0, 0, 0, 1; // SO3
Eigen::ColPivHouseholderQR<Eigen::MatrixXd> solver(compo_in_oxyde);
compo_species = solver.solve(compo_oxyde);
}
Event Timeline
Log In to Comment