chemmisol 0.1
Loading...
Searching...
No Matches
Classes | Public Member Functions | List of all members
chemmisol::ChemicalSystem Class Reference

#include <system.h>

Public Member Functions

 ChemicalSystem (const ChemicalSystem &other)
 
const ChemicalSystemoperator= (const ChemicalSystem &other)=delete
 
 ChemicalSystem (ChemicalSystem &&other)=delete
 
const ChemicalSystemoperator= (ChemicalSystem &&other)=delete
 
 ChemicalSystem ()=default
 
 ChemicalSystem (double solid_concentration, double specific_surface_area, double site_concentration)
 
void setUp ()
 
const ReactionaddReaction (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 ReactiongetReaction (const std::string &name) const
 
const ChemicalSpeciesgetSpecies (const std::string &name) const
 
const ChemicalSpeciesgetSpecies (const std::size_t &index) const
 
const ChemicalComponentgetComponent (const std::string &name) const
 
const ChemicalComponentgetComponent (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 ChemicalSpeciesReagentgetSpeciesReagent (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)
 

Detailed Description

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().

See also
ChemicalComponent
ChemicalSpecies
Reaction
Examples
basic_chemical_system/main.cpp.

Constructor & Destructor Documentation

◆ ChemicalSystem() [1/3]

chemmisol::ChemicalSystem::ChemicalSystem ( const ChemicalSystem other)

ChemicalSystem copy constructor.

◆ ChemicalSystem() [2/3]

chemmisol::ChemicalSystem::ChemicalSystem ( )
default

Defines a pure solution model, with only AQUEOUS components.

Trying to force an addComponent() call with MINERAL phase will yield undefined behaviors.

◆ ChemicalSystem() [3/3]

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.

Parameters
solid_concentrationQuantity of mineral in suspension in the solution, usually expressed in g/l.
specific_surface_areaSurface of the solid in contact with the solution per unit of mass, usually expressed in m2/g.
site_concentrationQuantity of sites per unit of surface in contact with the solution, usually expressed as entities/nm2.

Member Function Documentation

◆ setUp()

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.

Exceptions
EmptyReagentsif a reaction with an empty reagents list is processed.
MissingProducedSpeciesInReactionif a reaction is missing a produced species.
TooManyProducedSpeciesInReactionif several produced species seems to be specified and no unique produced species can be identified.

◆ addReaction()

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.

Parameters
nameReaction name.
log_KEquilibrium constant.
reagentsReagents of the reaction.
Returns
Reference to the new Reaction.
Examples
basic_chemical_system/main.cpp

◆ addComponent() [1/2]

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.

Parameters
nameName of the component.
total_concentrationInitial total concentration.
Examples
basic_chemical_system/main.cpp

◆ addComponent() [2/2]

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:

  • AQUEOUS: AqueousComponent
  • MINERAL: MineralComponent
Parameters
nameName of the component.
phaseChemical phase of the component.
total_concentrationInitial total concentration for AQUEOUS species, initial total molar fraction for MINERAL species.
Examples
basic_chemical_system/main.cpp
Exceptions
InvalidMineralSpeciesWithUndefinedSitesCountif 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.

◆ fixComponent() [1/2]

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.

Parameters
nameName of the component.
total_concentrationFixed total concentration.

◆ fixComponent() [2/2]

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:

  • AQUEOUS: AqueousComponent
  • MINERAL: MineralComponent The concentration of the Component is then fixed to the specified concentration.
Parameters
nameName of the component.
phaseChemical phase of the component.
total_concentrationFixed concentration for AQUEOUS species, fixed molar fraction for MINERAL species.
Exceptions
InvalidMineralSpeciesWithUndefinedSitesCountif 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.

◆ addSolvent()

void chemmisol::ChemicalSystem::addSolvent ( const std::string &  name)

Adds a new Solvent to the system.

Parameters
nameName of the solvent.
Examples
basic_chemical_system/main.cpp.

◆ initPH() [1/2]

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.

Parameters
pHInitial pH.
h_component_nameName of the component used to compute the pH.

◆ initPH() [2/2]

void chemmisol::ChemicalSystem::initPH ( double  pH)
inline

Initializes the pH in the default H+ component.

See also
initPH(double pH, const std::string& h_component_name)

◆ fixPH() [1/2]

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.

Parameters
pHFixed pH value.
h_component_nameName of the component used to compute the pH.
Examples
basic_chemical_system/main.cpp.

◆ fixPH() [2/2]

void chemmisol::ChemicalSystem::fixPH ( double  pH)
inline

Fixes the pH in the default H+ component.

See also
fixPH(double pH, const std::string& h_component_name)

◆ getPH() [1/2]

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.

◆ getPH() [2/2]

double chemmisol::ChemicalSystem::getPH ( ) const
inline

Returns the pH of the chemical system from the default H+ component.

◆ sitesQuantity()

double chemmisol::ChemicalSystem::sitesQuantity ( ) const
inline

Returns the quantity of mineral sites currently defined in the system. Might be null if the system is aqueous.

Returns
quantity of sites in mol

◆ setTotalQuantity()

void chemmisol::ChemicalSystem::setTotalQuantity ( const ChemicalComponent component,
double  quantity 
)

Sets the total quantity of the specified component.

Warning
The quantity of species produced from the specified component, including the species associated to the component, are not updated by this method. The solveEquilibrium() method must be run to update species quantities according the new equilibrium state and the mass conservation law.

◆ setTotalConcentration()

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.

Warning
The quantity of species produced from the specified component, including the species associated to the component, are not updated by this method. The solveEquilibrium() method must be run to update species quantities according the new equilibrium state and the mass conservation law.

◆ getReaction()

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().

Parameters
nameName of the reaction

◆ getSpecies() [1/3]

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.

Parameters
nameName of the component

◆ getSpecies() [2/3]

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.

Parameters
indexIndex of the component

◆ getComponent() [1/2]

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().

Parameters
nameName of the component.

◆ getComponent() [2/2]

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().

Parameters
indexIndex of the component.

◆ isComponent() [1/2]

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.

Parameters
species_nameName of the species to check.
Returns
True if the species is associated to a component.

◆ isComponent() [2/2]

bool chemmisol::ChemicalSystem::isComponent ( const ChemicalSpecies species) const

Checks if the species that corresponds to the specified name is associated to a component.

Parameters
speciesSpecies to check.
Returns
True if the species is associated to a component.

◆ getComponents()

const std::vector< std::unique_ptr< ChemicalComponent > > & chemmisol::ChemicalSystem::getComponents ( ) const
inline

Returns references to all the ChemicalComponents available in the system.

◆ getSpecies() [3/3]

const std::vector< std::unique_ptr< ChemicalSpecies > > & chemmisol::ChemicalSystem::getSpecies ( ) const
inline

Returns references to all the ChemicalSpecies available in the system.

◆ getReactions()

const std::vector< std::unique_ptr< Reaction > > & chemmisol::ChemicalSystem::getReactions ( ) const
inline

Returns references to all the Reactions added to the system.

◆ getComponentReagents()

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.

Parameters
reactionReference to a reaction added to the system.
Returns
List of reagents bound to ChemicalComponents of this chemical system.

◆ getSpeciesReagent()

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.

Parameters
reactionReference to a reaction added to the system.
Returns
Produced species of the reaction, bound to a ChemicalSpecies of this chemical system.

◆ getReactionMatrix()

const std::vector< std::vector< double > > & chemmisol::ChemicalSystem::getReactionMatrix ( ) const
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.

Example

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().

Returns
Matrix representing the reactions of the chemical system.

◆ proceed()

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).

Parameters
reactionReaction to proceed.
extentRequired extent of the the reaction (positive or negative).

◆ distanceToEquilibrium() [1/2]

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.

Parameters
activitiesActivities 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.
reactionReaction in this chemical system.
Returns
Distance to the equilibrium state of the specified reaction.

◆ degree()

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.

◆ distanceToEquilibrium() [2/2]

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&)

Parameters
reactionReaction in this chemical system.
Returns
Distance to the equilibrium state of the specified reaction.

◆ solveEquilibrium() [1/2]

void chemmisol::ChemicalSystem::solveEquilibrium ( )

Solves the equilibrium of the system using the default_solver.

Concentrations of all components are updated accordingly upon return.

Exceptions
EmptyReagentsif a reaction with an empty reagents list is processed.
MissingProducedSpeciesInReactionif a reaction is missing a produced species.
TooManyProducedSpeciesInReactionif several produced species seems to be specified and no unique produced species can be identified.
Examples
basic_chemical_system/main.cpp.

◆ solveEquilibrium() [2/2]

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.

Exceptions
EmptyReagentsif a reaction with an empty reagents list is processed.
MissingProducedSpeciesInReactionif a reaction is missing a produced species.
TooManyProducedSpeciesInReactionif several produced species seems to be specified and no unique produced species can be identified.

◆ reactionQuotient() [1/2]

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.

Parameters
reactionReaction in this chemical system.
Returns
Reaction quotient of the reaction.

◆ reactionQuotient() [2/2]

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().

See also
reactionQuotient(const Reaction&)

◆ massConservationLaw() [1/2]

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.

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).

Note

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.

Parameters
activitiesActivities of all species in the system, so that activities[i] corresponds to the activity of the species with index i.
resultVector 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.

◆ massConservationLaw() [2/2]

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.

Parameters
resultVector 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.

◆ getMaxIteration()

std::size_t chemmisol::ChemicalSystem::getMaxIteration ( ) const
inline

Returns the maximum count of iterations allowed for the equilibrium solver (200 by default).

◆ setMaxIteration()

void chemmisol::ChemicalSystem::setMaxIteration ( const std::size_t &  max_iteration)
inline

Sets the maximum count of iterations allowed for the equilibrium solver (200 by default).


The documentation for this class was generated from the following file: