chemmisol 0.1
|
#include <system.h>
Public Member Functions | |
ChemicalSystem (const ChemicalSystem &other) | |
const ChemicalSystem & | operator= (const ChemicalSystem &other)=delete |
ChemicalSystem (ChemicalSystem &&other)=delete | |
const ChemicalSystem & | operator= (ChemicalSystem &&other)=delete |
ChemicalSystem ()=default | |
ChemicalSystem (double solid_concentration, double specific_surface_area, double site_concentration) | |
void | setUp () |
const Reaction & | addReaction (const std::string &name, double log_K, const std::vector< Reagent > &reagents) |
void | addComponent (const std::string &name, double total_concentration) |
void | addComponent (const std::string &name, Phase phase, double total_concentration) |
void | fixComponent (const std::string &name, double total_concentration) |
void | fixComponent (const std::string &name, Phase phase, double total_concentration) |
void | addSolvent (const std::string &name) |
void | initPH (double pH, const std::string &h_component_name) |
void | initPH (double pH) |
void | fixPH (double pH, const std::string &h_component_name) |
void | fixPH (double pH) |
double | getPH (const std::string &h_component_name) const |
double | getPH () const |
double | sitesQuantity () const |
void | setTotalQuantity (const ChemicalComponent &component, double quantity) |
void | setTotalConcentration (const ChemicalComponent &component, double concentration) |
const Reaction & | getReaction (const std::string &name) const |
const ChemicalSpecies & | getSpecies (const std::string &name) const |
const ChemicalSpecies & | getSpecies (const std::size_t &index) const |
const ChemicalComponent & | getComponent (const std::string &name) const |
const ChemicalComponent & | getComponent (const std::size_t &index) const |
bool | isComponent (const std::string &species_name) const |
bool | isComponent (const ChemicalSpecies &species) const |
const std::vector< std::unique_ptr< ChemicalComponent > > & | getComponents () const |
const std::vector< std::unique_ptr< ChemicalSpecies > > & | getSpecies () const |
const std::vector< std::unique_ptr< Reaction > > & | getReactions () const |
const std::vector< ComponentReagent > & | getComponentReagents (const Reaction &reaction) const |
const ChemicalSpeciesReagent & | getSpeciesReagent (const Reaction &reaction) const |
const std::vector< std::vector< double > > & | getReactionMatrix () const |
void | proceed (const Reaction &reaction, double extent) |
double | distanceToEquilibrium (const std::vector< double > &activities, const Reaction &reaction) const |
std::complex< double > | distanceToEquilibrium (const std::vector< std::complex< double > > &activities, const Reaction &reaction) const |
double | degree (const Reaction &reaction) const |
double | distanceToEquilibrium (const Reaction &reaction) const |
void | solveEquilibrium () |
void | solveEquilibrium (const solver::Solver &solver) |
double | reactionQuotient (const Reaction &reaction) const |
double | reactionQuotient (const std::string &name) const |
void | massConservationLaw (const std::vector< double > &activities, std::vector< double > &result) const |
void | massConservationLaw (const std::vector< std::complex< double > > &activities, std::vector< std::complex< double > > &result) const |
void | massConservationLaw (std::vector< double > &result) const |
std::size_t | getMaxIteration () const |
void | setMaxIteration (const std::size_t &max_iteration) |
A ChemicalSystem is defined as a set of Components that interact according to Reactions in order to produce ChemicalSpecies.
A ChemicalSystem can be used to define both a pure solution system (where all species are aqueous) as well as a mineral adsorption model.
Once components and reactions are specified, the solveEquilibrium() method can be used to run a solver and set the activities of all chemical species so that the law of conservation of mass and chemical equilibriums are satisfied.
All quantities must be specified using the chemmisol unit system, or directly in core units. See UNITS().
chemmisol::ChemicalSystem::ChemicalSystem | ( | const ChemicalSystem & | other | ) |
ChemicalSystem copy constructor.
|
default |
Defines a pure solution model, with only AQUEOUS components.
Trying to force an addComponent() call with MINERAL phase will yield undefined behaviors.
chemmisol::ChemicalSystem::ChemicalSystem | ( | double | solid_concentration, |
double | specific_surface_area, | ||
double | site_concentration | ||
) |
Defines an adsorption model that can handle both AQUEOUS and MINERAL species. See MineralComponent for more detailed information about how mineral parameters are defined and used.
The surface complex corresponds to the "free sites" species. It is automatically added as MineralComponent with a total molar fraction of 1.
solid_concentration | Quantity of mineral in suspension in the solution, usually expressed in g/l. |
specific_surface_area | Surface of the solid in contact with the solution per unit of mass, usually expressed in m2/g. |
site_concentration | Quantity of sites per unit of surface in contact with the solution, usually expressed as entities/nm2. |
void chemmisol::ChemicalSystem::setUp | ( | ) |
Sets up the system so that it is ready to solve.
The method is automatically called by solveEquilibrium(), but can be called by the user to inspect the state of the system before solving.
Reagents of reactions that have not been added with addSpecies(), addComponent() or fixComponent() are automatically added as chemical species with an initially null concentration.
EmptyReagents | if a reaction with an empty reagents list is processed. |
MissingProducedSpeciesInReaction | if a reaction is missing a produced species. |
TooManyProducedSpeciesInReaction | if several produced species seems to be specified and no unique produced species can be identified. |
const Reaction & chemmisol::ChemicalSystem::addReaction | ( | const std::string & | name, |
double | log_K, | ||
const std::vector< Reagent > & | reagents | ||
) |
Adds a new Reaction to the system.
See Reaction for the precise meaning of each parameter.
name | Reaction name. |
log_K | Equilibrium constant. |
reagents | Reagents of the reaction. |
void chemmisol::ChemicalSystem::addComponent | ( | const std::string & | name, |
double | total_concentration | ||
) |
Adds a new AqueousComponent to the chemical system, assuming that the Component is AQUEOUS by default.
name | Name of the component. |
total_concentration | Initial total concentration. |
void chemmisol::ChemicalSystem::addComponent | ( | const std::string & | name, |
Phase | phase, | ||
double | total_concentration | ||
) |
Adds a new Component to the chemical system, depending on the provided phase, the following Component is instantiated:
name | Name of the component. |
phase | Chemical phase of the component. |
total_concentration | Initial total concentration for AQUEOUS species, initial total molar fraction for MINERAL species. |
InvalidMineralSpeciesWithUndefinedSitesCount | if the phase is MINERAL but the sites quantity of the system is null. See the ChemicalSystem(double, double, double, const std::string&) to ensure the sites quantity is properly defined. |
void chemmisol::ChemicalSystem::fixComponent | ( | const std::string & | name, |
double | total_concentration | ||
) |
Adds a new AqueousComponent to the chemical system, assuming that the Component is AQUEOUS by default, and fixes it to the specified concentration.
name | Name of the component. |
total_concentration | Fixed total concentration. |
void chemmisol::ChemicalSystem::fixComponent | ( | const std::string & | name, |
Phase | phase, | ||
double | total_concentration | ||
) |
Adds a new Component to the chemical system, depending on the provided phase, the following Component is instantiated:
name | Name of the component. |
phase | Chemical phase of the component. |
total_concentration | Fixed concentration for AQUEOUS species, fixed molar fraction for MINERAL species. |
InvalidMineralSpeciesWithUndefinedSitesCount | if the phase is MINERAL but the sites quantity of the system is null. See the ChemicalSystem(double, double, double, const std::string&) to ensure the sites quantity is properly defined. |
void chemmisol::ChemicalSystem::addSolvent | ( | const std::string & | name | ) |
Adds a new Solvent to the system.
name | Name of the solvent. |
void chemmisol::ChemicalSystem::initPH | ( | double | pH, |
const std::string & | h_component_name | ||
) |
Initializes the pH of the chemical system, adding a component with the specified name and a concentration of 10^-pH.
Notice that this only consists in a start value, that is very likely to vary when the equilibrium is solved, except if fixPH() was called.
pH | Initial pH. |
h_component_name | Name of the component used to compute the pH. |
|
inline |
Initializes the pH in the default H+ component.
void chemmisol::ChemicalSystem::fixPH | ( | double | pH, |
const std::string & | h_component_name | ||
) |
Fixed the pH of the system to the specified value in the component with the specified name, that is created if it does not exist yet.
pH | Fixed pH value. |
h_component_name | Name of the component used to compute the pH. |
|
inline |
Fixes the pH in the default H+ component.
double chemmisol::ChemicalSystem::getPH | ( | const std::string & | h_component_name | ) | const |
Returns the pH of the chemical system from the component with the specified name.
|
inline |
Returns the pH of the chemical system from the default H+ component.
|
inline |
Returns the quantity of mineral sites currently defined in the system. Might be null if the system is aqueous.
void chemmisol::ChemicalSystem::setTotalQuantity | ( | const ChemicalComponent & | component, |
double | quantity | ||
) |
Sets the total quantity of the specified component.
void chemmisol::ChemicalSystem::setTotalConcentration | ( | const ChemicalComponent & | component, |
double | concentration | ||
) |
Sets the total quantity of the specified component to a quantity that corresponds to the quantity of its associated species at the provided concentration.
const Reaction & chemmisol::ChemicalSystem::getReaction | ( | const std::string & | name | ) | const |
Gets the reaction named name
.
The behavior of the method is unspecified if the name does not correspond to any reaction added with addReaction().
name | Name of the reaction |
const ChemicalSpecies & chemmisol::ChemicalSystem::getSpecies | ( | const std::string & | name | ) | const |
Gets the chemical species named name
.
The behavior of the method is unspecified if the name does not correspond to any chemical species in the system. Valid species include species added explicitly by addComponent(), fixComponent() or addSolvent(), and species added implicitly as reagents of reactions by the setUp() method.
name | Name of the component |
const ChemicalSpecies & chemmisol::ChemicalSystem::getSpecies | ( | const std::size_t & | index | ) | const |
Gets the component with the specified index, that can be retrieved from an existing component with Component::getIndex().
The behavior of the method is unspecified if the index does not correspond to any chemical species in the system. Valid species include species added explicitly by addComponent(), fixComponent() or addSolvent(), and species added implicitly as reagents of reactions by the setUp() method.
index | Index of the component |
const ChemicalComponent & chemmisol::ChemicalSystem::getComponent | ( | const std::string & | name | ) | const |
Gets the component named name
.
The behavior of the method is unspecified if the name does not correspond to any component in the system. Valid components include components added explicitly with addComponent(), fixComponent() or addSolvent().
name | Name of the component. |
const ChemicalComponent & chemmisol::ChemicalSystem::getComponent | ( | const std::size_t & | index | ) | const |
Gets the component with the specified index, that can be retrieved from an existing component with Component::getIndex().
The behavior of the method is unspecified if the index does not correspond to any component in the system. Valid components include components added explicitly with addComponent(), fixComponent() or addSolvent().
index | Index of the component. |
bool chemmisol::ChemicalSystem::isComponent | ( | const std::string & | species_name | ) | const |
Checks if the species that corresponds to the specified name is associated to a component.
The behavior of the method is unspecified if the name does not correspond to any chemical species in the system. Valid species include species added explicitly by addComponent(), fixComponent() or addSolvent(), and species added implicitly as reagents of reactions by the setUp() method.
species_name | Name of the species to check. |
bool chemmisol::ChemicalSystem::isComponent | ( | const ChemicalSpecies & | species | ) | const |
Checks if the species that corresponds to the specified name is associated to a component.
species | Species to check. |
|
inline |
Returns references to all the ChemicalComponents available in the system.
|
inline |
Returns references to all the ChemicalSpecies available in the system.
|
inline |
Returns references to all the Reactions added to the system.
const std::vector< ComponentReagent > & chemmisol::ChemicalSystem::getComponentReagents | ( | const Reaction & | reaction | ) | const |
Gets the list of compiled component reagents of the specified reaction.
The behavior of the method is unspecified if setUp() has not been called since the specified reaction was added to the system.
reaction | Reference to a reaction added to the system. |
const ChemicalSpeciesReagent & chemmisol::ChemicalSystem::getSpeciesReagent | ( | const Reaction & | reaction | ) | const |
Gets the compiled produced species of the specified reaction.
The behavior of the method is unspecified if setUp() has not been called since the specified reaction was added to the system.
reaction | Reference to a reaction added to the system. |
|
inline |
Returns the reaction matrix describing the system, that should have been previously computed by the setUp() method.
Each row corresponds to a Reaction of the system, each column corresponds to a ChemicalSpecies of the system, and the coefficient at (i, j) correspond to the stoichiometric coefficient of the ChemicalSpecies j in the Reaction i. Notice that the coefficients might be 0.
Considering the following chemical system:
HO- : H2O <-> H+ + HO- NaCl: Na+ + Cl- <-> NaCl NaOH: Na+ + H2O <-> NaOH + H+
The reaction matrix can be defined as follows:
H2O | HO- | H+ | NaCl | Na+ | Cl- | NaOH | |
---|---|---|---|---|---|---|---|
HO- | 1 | -1 | -1 | 0 | 0 | 0 | 0 |
NaCl | 0 | 0 | 0 | -1 | 1 | 1 | 0 |
NaOH | 1 | 0 | -1 | 0 | 1 | 0 | -1 |
The actual indexes of each reaction and component correspond to Reaction::getIndex() and ChemicalSpecies::getIndex().
void chemmisol::ChemicalSystem::proceed | ( | const Reaction & | reaction, |
double | extent | ||
) |
Updates the quantities of species in the system resulting from the specified extent of the reaction.
It is the responsibility of the user to ensure that extent is in a valid range, according to limiting factors of the reaction.
According to the Chemmisol convention about stoichiometric coefficients (see Reaction), a positive extent increases the quantity of reactants, and a negative extent increases the quantity of products (since products are associated to negative coefficients, by convention).
reaction | Reaction to proceed. |
extent | Required extent of the the reaction (positive or negative). |
double chemmisol::ChemicalSystem::distanceToEquilibrium | ( | const std::vector< double > & | activities, |
const Reaction & | reaction | ||
) | const |
Computes the distance between the system state represented by the specified activities and the equilibrium state of the specified reaction.
Considering the following reaction equation (see Reaction for more details) with an equilibirum constant K:
n A + m B - k C - l D <-> 0
the distance to the equilibrium is computed as:
d = K * a(C)^k * a(D)^l - (a(A)^n + a(B)^m)
where a(S) denotes the activity of the species S, with a(S)^n set to 0 if a(S)=0.
The distance to the equilibrium thus constitutes a polynomial equation where variables represent the activity of each species, and the equilibrium of the system constitutes roots of the polynom. The total degree of this polynom can be queried with the degree() method. In this case, the total degree would be max(k+l, n+m)
.
Such distance might be positive or negative, depending on which side of the equilibrium the current state is.
activities | Activities considered to compute the distance to equilibrium. Activities must be specified so that activities[species ->getIndex()] is equal to the activity of species , for any species corresponding to a reagent of the reaction. |
reaction | Reaction in this chemical system. |
double chemmisol::ChemicalSystem::degree | ( | const Reaction & | reaction | ) | const |
Returns the total degree of the polynomial equation representing the "distance to equilibrium equation".
See distanceToEquilibrium() for an example.
double chemmisol::ChemicalSystem::distanceToEquilibrium | ( | const Reaction & | reaction | ) | const |
Computes the distance between the current state of the system and the equilibrium state of the specified reaction.
See distanceToEquilibrium(const std::vector<double>&, const Reaction&)
reaction | Reaction in this chemical system. |
void chemmisol::ChemicalSystem::solveEquilibrium | ( | ) |
Solves the equilibrium of the system using the default_solver.
Concentrations of all components are updated accordingly upon return.
EmptyReagents | if a reaction with an empty reagents list is processed. |
MissingProducedSpeciesInReaction | if a reaction is missing a produced species. |
TooManyProducedSpeciesInReaction | if several produced species seems to be specified and no unique produced species can be identified. |
void chemmisol::ChemicalSystem::solveEquilibrium | ( | const solver::Solver & | solver | ) |
Solves the equilibrium of the system using the provided solver.
Concentrations of all components are updated accordingly upon return.
EmptyReagents | if a reaction with an empty reagents list is processed. |
MissingProducedSpeciesInReaction | if a reaction is missing a produced species. |
TooManyProducedSpeciesInReaction | if several produced species seems to be specified and no unique produced species can be identified. |
double chemmisol::ChemicalSystem::reactionQuotient | ( | const Reaction & | reaction | ) | const |
Computes the reaction quotient of the specified reaction.
By convention, the products of the reaction (i.e. reagents with a negative coefficient) form the numerous of the quotient.
Returns a NaN value if the activity of one of the reactants is null.
reaction | Reaction in this chemical system. |
double chemmisol::ChemicalSystem::reactionQuotient | ( | const std::string & | name | ) | const |
Computes the reaction quotient of the reaction corresponding to the specified name.
The behavior of the method is unspecified if the name does not correspond to any reaction in the system. Valid reactions include reactions added explicitly with addReaction().
void chemmisol::ChemicalSystem::massConservationLaw | ( | const std::vector< double > & | activities, |
std::vector< double > & | result | ||
) | const |
Computes the value of the law of conservation of mass for a system state corresponding to the provided species activity.
The actual total quantity of a component is defined as the sum of the quantity of species associated to the component itself and the quantities of all species compound from this component, weighted by the stoichiometric coefficient associated to the component in the reaction producing exactly one unit of the compound. See the ChemicalComponent documentation for a detailed example.
Let's define the same system as the example of the ChemicalSystem documentation, with the following components: H+ and PO4.
The following reactions can be defined to introduce the OH- and H3PO4 produced species as follows:
H2O <-> OH- + H+ PO4 + 3H+ <-> H3PO4
The actual (or "output") total quantity N of the PO4 and H+ components are then defined as follows:
N(PO4) = n(PO4) + n(H3PO4) N(H+) = n(H+) - n(OH-) + 3 * n(H3PO4)
Then, let's define Q(C) as the input total quantity of the component C, that can be set with ChemicalComponent::setTotalQuantity(). The law of conservation of mass states that at any time in the system, the following relations should hold:
Q(PO4) = N(PO4) Q(H+) = N(H+)
and so:
n(PO4) + n(H3PO4) - Q(PO4) = 0 n(H+) - n(OH-) + 3 * n(H3PO4) - Q(H+) = 0
Those constraints are thus added to the equation system solved by the solveEquilibrium() method, next to the chemical equilibrium constraints for each reaction.
The massConservationLaw() method then returns a vector such that result[i] = N(C[i])-Q(C[i])
where C[i]
represents the component for index i, except if the component C[i]
is fixed. In this case, result[i] = 0 (see the note below).
The usage of the mass conservation law is only relevant for components that are not fixed. Indeed, the notion of "fixed
component" cannot exist in a closed real system, since the extent of any reaction will modify the quantity of fixed species, including solvents. In order to apply the law of conservation of mass to a fixed component, it is necessary to implicitly define a "buffer" (for example, a pH buffer solution) that can consume and produce any desired quantity of the species associated to the fixed component to keep the quantity of the fixed species to the desired value . The law of mass conservation is then rather expressed as a mass balance equation of the form Input + Generation = Output + Accumulation + Consumption
. For example, the mass conservation of H2O in the above system should then be expressed as:
N(H2O) + H2O production = n(H2O) + n(OH-) + H2O consumption
so that H2O production and H2O consumption are set to any value that guarantees that n(H2O) is fixed to the desired value, whatever the quantity of OH- (and so the extent of the H20 <-> H+ + OH- reaction) is.
Rather than trying to simulate the action of a buffer, using a fixed component implies the existence of a perfect buffer, what is mathematically represented by the fact that in any case Q(C)-N(C)=0
for any fixed species. In consequence, no constraint is added to the solved system for fixed species and result[index]=0
for any index corresponding to a fixed component.
activities | Activities of all species in the system, so that activities[i] corresponds to the activity of the species with index i. |
result | Vector in which massConservationLaw results for all components in the system are written, so that result[i] corresponds to the result of the mass conservation law for the component with index i. |
void chemmisol::ChemicalSystem::massConservationLaw | ( | std::vector< double > & | result | ) | const |
Computes the value of the law of conservation of mass for the current system state.
See massConservationLaw(activities, result) for detailed information.
result | Vector in which massConservationLaw results for all components in the system are written, so that result[i] corresponds to the result of the mass conservation law for the component with index i. |
|
inline |
Returns the maximum count of iterations allowed for the equilibrium solver (200 by default).
|
inline |
Sets the maximum count of iterations allowed for the equilibrium solver (200 by default).