Class ChemEquil implements a chemical equilibrium solver for single-phase solutions. More...

#include <ChemEquil.h>

Collaboration diagram for ChemEquil:

Public Member Functions
	ChemEquil ()
	Default Constructor.

	ChemEquil (thermo_t &s)
	Constructor combined with the initialization function.

int	equilibrate (thermo_t &s, const char *XY, bool useThermoPhaseElementPotentials=false, int loglevel=0)
	Equilibrate a phase, holding the elemental composition fixed at the initial value found within the ThermoPhase object.

int	equilibrate (thermo_t &s, const char *XY, vector_fp &elMoles, bool useThermoPhaseElementPotentials=false, int loglevel=0)
	Compute the equilibrium composition for 2 specified properties and the specified element moles.

const vector_fp &	elementPotentials () const

Public Attributes
EquilOpt	options
	Options controlling how the calculation is carried out.

Protected Member Functions
doublereal	nAtoms (size_t k, size_t m) const
	number of atoms of element m in species k.

void	initialize (thermo_t &s)
	Prepare for equilibrium calculations.

void	setToEquilState (thermo_t &s, const vector_fp &x, doublereal t)
	Set mixture to an equilibrium state consistent with specified element potentials and temperature.

int	setInitialMoles (thermo_t &s, vector_fp &elMoleGoal, int loglevel=0)
	Estimate the initial mole numbers.

int	estimateElementPotentials (thermo_t &s, vector_fp &lambda, vector_fp &elMolesGoal, int loglevel=0)
	Generate a starting estimate for the element potentials.

int	estimateEP_Brinkley (thermo_t &s, vector_fp &lambda, vector_fp &elMoles)
	Do a calculation of the element potentials using the Brinkley method, p.

int	dampStep (thermo_t &s, vector_fp &oldx, double oldf, vector_fp &grad, vector_fp &step, vector_fp &x, double &f, vector_fp &elmols, double xval, double yval)

void	equilResidual (thermo_t &s, const vector_fp &x, const vector_fp &elmtotal, vector_fp &resid, double xval, double yval, int loglevel=0)
	Evaluates the residual vector F, of length mm.

void	equilJacobian (thermo_t &s, vector_fp &x, const vector_fp &elmols, DenseMatrix &jac, double xval, double yval, int loglevel=0)

void	adjustEloc (thermo_t &s, vector_fp &elMolesGoal)

void	update (const thermo_t &s)
	update internally stored state information.

double	calcEmoles (thermo_t &s, vector_fp &x, const double &n_t, const vector_fp &Xmol_i_calc, vector_fp &eMolesCalc, vector_fp &n_i_calc, double pressureConst)
	Given a vector of dimensionless element abundances, this routine calculates the moles of the elements and the moles of the species.

Protected Attributes
thermo_t *	m_phase
	Pointer to the ThermoPhase object used to initialize this object.

size_t	m_mm

size_t	m_kk

size_t	m_skip

size_t	m_nComponents
	This is equal to the rank of the stoichiometric coefficient matrix when it is computed.

std::auto_ptr < PropertyCalculator< thermo_t > >	m_p1

std::auto_ptr < PropertyCalculator< thermo_t > >	m_p2

vector_fp	m_molefractions
	Current value of the mole fractions in the single phase.

vector_fp	m_lambda
	Current value of the dimensional element potentials -> length = m_mm.

doublereal	m_elementTotalSum

vector_fp	m_elementmolefracs

vector_fp	m_reswork

vector_fp	m_jwork1

vector_fp	m_jwork2

vector_fp	m_comp

doublereal	m_temp

doublereal	m_dens

doublereal	m_p0

size_t	m_eloc
	Index of the element id corresponding to the electric charge of each species.

vector_fp	m_startSoln

vector_fp	m_grt

vector_fp	m_mu_RT

vector_fp	m_muSS_RT
	Dimensionless values of the gibbs free energy for the standard state of each species, at the temperature and pressure of the solution (the star standard state).

std::vector< size_t >	m_component

double	m_elemFracCutoff

bool	m_doResPerturb

std::vector< size_t >	m_orderVectorElements

std::vector< size_t >	m_orderVectorSpecies

Detailed Description

Class ChemEquil implements a chemical equilibrium solver for single-phase solutions.

It is a "non-stoichiometric" solver in the terminology of Smith and Missen, meaning that every intermediate state is a valid chemical equilibrium state, but does not necessarily satisfy the element constraints. In contrast, the solver implemented in class MultiPhaseEquil uses a "stoichiometric" algorithm, in which each intermediate state satisfies the element constraints but is not a state of chemical equilibrium. Non-stoichiometric methods are faster when they converge, but stoichiometric ones tend to be more robust and can be used also for problems with multiple condensed phases. As expected, the ChemEquil solver is faster than MultiPhaseEquil for many single-phase equilibrium problems (particularly if there are only a few elements but very many species), but can be less stable. Problem situations include low temperatures where only a few species have non-zero mole fractions, precisely stoichiometric compositions (e.g. 2 H2 + O2). In general, if speed is important, this solver should be tried first, and if it fails then use MultiPhaseEquil.

Definition at line 98 of file ChemEquil.h.

Constructor & Destructor Documentation

ChemEquil ( )

Default Constructor.

Definition at line 67 of file ChemEquil.cpp.

ChemEquil ( thermo_t & s )

Constructor combined with the initialization function.

This constructor initializes the ChemEquil object with everything it needs to start solving equilibrium problems.

Parameters

s	ThermoPhase object that will be used in the equilibrium calls.

Definition at line 79 of file ChemEquil.cpp.

References ChemEquil::initialize().

Member Function Documentation

int equilibrate	(	thermo_t &	s,
		const char *	XY,
		bool	useThermoPhaseElementPotentials = `false`,
		int	loglevel = `0`
	)

Equilibrate a phase, holding the elemental composition fixed at the initial value found within the ThermoPhase object.

The value of 2 specified properties are obtained by querying the ThermoPhase object. The properties must be already contained within the current thermodynamic state of the system.

Definition at line 452 of file ChemEquil.cpp.

References ChemEquil::initialize(), Phase::nElements(), and ChemEquil::update().

Referenced by Cantera::equilibrate(), and Cantera::vcs_equilibrate().

int equilibrate	(	thermo_t &	s,
		const char *	XYstr,
		vector_fp &	elMolesGoal,
		bool	useThermoPhaseElementPotentials = `false`,
		int	loglevel = `0`
	)

Compute the equilibrium composition for 2 specified properties and the specified element moles.

elMoles = specified vector of element abundances.

The 2 specified properties are obtained by querying the ThermoPhase object. The properties must be already contained within the current thermodynamic state of the system.

Return variable: Successful returns are indicated by a return value of 0. Unsuccessful returns are indicated by a return value of -1 for lack of convergence or -3 for a singular jacobian.

Definition at line 480 of file ChemEquil.cpp.

References Cantera::_equilflag(), EquilOpt::absElemTol, Cantera::addLogEntry(), Cantera::beginLogGroup(), DATA_PTR, Cantera::dot(), Phase::elementName(), Phase::elementNames(), Cantera::endLogGroup(), ChemEquil::equilResidual(), ChemEquil::estimateElementPotentials(), ChemEquil::estimateEP_Brinkley(), Cantera::fp2str(), Cantera::GasConstant, ThermoPhase::getElementPotentials(), ChemEquil::initialize(), Cantera::int2str(), EquilOpt::iterations, DenseMatrix::leftMult(), ChemEquil::m_eloc, ChemEquil::m_lambda, ChemEquil::m_nComponents, ckr::max(), EquilOpt::maxIterations, ThermoPhase::maxTemp(), ckr::min(), ThermoPhase::minTemp(), Phase::moleFraction(), Phase::nElements(), Cantera::npos, Phase::nSpecies(), ChemEquil::options, EquilOpt::relTolerance, Phase::restoreState(), CanteraError::save(), Phase::saveState(), Cantera::scale(), ThermoPhase::setElementPotentials(), ChemEquil::setInitialMoles(), Phase::setMoleFractions(), Phase::setTemperature(), ChemEquil::setToEquilState(), Cantera::solve(), Phase::temperature(), ChemEquil::update(), Cantera::writelog(), and Cantera::writelogf().

doublereal nAtoms	(	size_t	k,
		size_t	m
	)		const

inlineprotected

number of atoms of element m in species k.

Definition at line 143 of file ChemEquil.h.

Referenced by ChemEquil::calcEmoles(), ChemEquil::estimateElementPotentials(), ChemEquil::estimateEP_Brinkley(), ChemEquil::setToEquilState(), and ChemEquil::update().

void initialize ( thermo_t & s )

protected

Prepare for equilibrium calculations.

Parameters

s	object representing the solution phase.

Definition at line 97 of file ChemEquil.cpp.

References Phase::atomicWeight(), Phase::elementName(), Cantera::fp2str(), ChemEquil::m_eloc, ChemEquil::m_lambda, ChemEquil::m_molefractions, ChemEquil::m_muSS_RT, ChemEquil::m_nComponents, ChemEquil::m_phase, Phase::nAtoms(), Phase::nElements(), Cantera::npos, Phase::nSpecies(), ThermoPhase::refPressure(), Phase::speciesName(), and Cantera::writelog().

Referenced by ChemEquil::ChemEquil(), and ChemEquil::equilibrate().

void setToEquilState	(	thermo_t &	s,
		const vector_fp &	lambda_RT,
		doublereal	t
	)

protected

Set mixture to an equilibrium state consistent with specified element potentials and temperature.

Parameters

lambda_RT	vector of non-dimensional element potentials \[ \lambda_m/RT \] .
t	temperature in K.

Definition at line 182 of file ChemEquil.cpp.

References DATA_PTR, ChemEquil::nAtoms(), Phase::setTemperature(), ThermoPhase::setToEquilState(), and ChemEquil::update().

Referenced by ChemEquil::equilibrate(), and ChemEquil::equilResidual().

int setInitialMoles	(	thermo_t &	s,
		vector_fp &	elMoleGoal,
		int	loglevel = `0`
	)

protected

Estimate the initial mole numbers.

This version borrows from the MultiPhaseEquil solver.

Definition at line 238 of file ChemEquil.cpp.

References Cantera::addLogEntry(), MultiPhase::addPhase(), Cantera::beginLogGroup(), PrintCtrl::cropAbs10(), PrintCtrl::CT_OFF_GLOBALOBEY, Phase::elementName(), Cantera::endLogGroup(), MultiPhase::init(), ChemEquil::m_nComponents, Phase::moleFraction(), ThermoPhase::pressure(), CanteraError::save(), Phase::speciesName(), Phase::temperature(), ChemEquil::update(), Cantera::writelog(), and Cantera::writelogf().

Referenced by ChemEquil::equilibrate().

int estimateElementPotentials	(	thermo_t &	s,
		vector_fp &	lambda,
		vector_fp &	elMolesGoal,
		int	loglevel = `0`
	)

protected

Generate a starting estimate for the element potentials.

Definition at line 309 of file ChemEquil.cpp.

References Cantera::addLogEntry(), MultiPhase::addPhase(), Cantera::BasisOptimize(), Cantera::beginLogGroup(), PrintCtrl::cropAbs10(), PrintCtrl::CT_OFF_GLOBALOBEY, DATA_PTR, Phase::elementName(), Cantera::ElemRearrange(), Cantera::endLogGroup(), Cantera::GasConstant, ThermoPhase::getChemPotentials(), Phase::getMoleFractions(), MultiPhase::init(), ChemEquil::m_nComponents, ChemEquil::nAtoms(), Phase::nSpecies(), ThermoPhase::pressure(), Cantera::scale(), Phase::setMoleFractions(), Cantera::solve(), Phase::speciesName(), Phase::temperature(), Cantera::writelog(), and Cantera::writelogf().

Referenced by ChemEquil::equilibrate().

int estimateEP_Brinkley	(	thermo_t &	s,
		vector_fp &	x,
		vector_fp &	elMoles
	)

protected

Do a calculation of the element potentials using the Brinkley method, p.

129 Smith and Missen.

We have found that the previous estimate may not be good enough to avoid drastic numerical issues associated with the use of a numerically generated jacobian used in the main algorithm.

The Brinkley algorithm, here, assumes a constant T, P system and uses a linearized analytical Jacobian that turns out to be very stable even given bad initial guesses.

The pressure and temperature to be used are in the ThermoPhase object input into the routine.

The initial guess for the element potentials used by this routine is taken from the input vector, x.

elMoles is the input element abundance vector to be matched.

Nonideal phases are handled in principle. This is done by calculating the activity coefficients and adding them into the formula in the correct position. However, these are treated as a rhs contribution only. Therefore, convergence might be a problem. This has not been tested. Also molality based unit systems aren't handled.

On return, int return value contains the success code: 0 - successful 1 - unsuccessful, max num iterations exceeded -3 - unsuccessful, singular jacobian

NOTE: update for activity coefficients.

Definition at line 1326 of file ChemEquil.cpp.

References EquilOpt::absElemTol, Cantera::addLogEntry(), ChemEquil::calcEmoles(), DATA_PTR, Phase::elementName(), Phase::elementNames(), ThermoPhase::getActivityCoefficients(), ThermoPhase::getGibbs_RT(), Phase::getMoleFractions(), ChemEquil::m_eloc, ChemEquil::m_muSS_RT, ChemEquil::m_nComponents, EquilOpt::maxIterations, ChemEquil::nAtoms(), Cantera::npos, ChemEquil::options, ThermoPhase::pressure(), EquilOpt::relTolerance, Phase::restoreState(), CanteraError::save(), Phase::saveState(), Phase::setMoleFractions(), ThermoPhase::setPressure(), Cantera::solve(), Phase::speciesName(), Phase::temperature(), Cantera::writelog(), and Cantera::writelogf().

Referenced by ChemEquil::equilibrate().

void equilResidual	(	thermo_t &	s,
		const vector_fp &	x,
		const vector_fp &	elmtotal,
		vector_fp &	resid,
		double	xval,
		double	yval,
		int	loglevel = `0`
	)

protected

Evaluates the residual vector F, of length mm.

Definition at line 1120 of file ChemEquil.cpp.

References Cantera::addLogEntry(), Cantera::beginLogGroup(), PrintCtrl::cropAbs10(), PrintCtrl::CT_OFF_GLOBALOBEY, Phase::elementName(), Cantera::endLogGroup(), Cantera::fp2str(), ChemEquil::m_eloc, ChemEquil::m_nComponents, ChemEquil::setToEquilState(), Cantera::writelog(), and Cantera::writelogf().

Referenced by ChemEquil::equilibrate().

void update ( const thermo_t & s )

protected

update internally stored state information.

Definition at line 206 of file ChemEquil.cpp.

References DATA_PTR, Phase::density(), Cantera::fp2str(), Phase::getMoleFractions(), ChemEquil::m_molefractions, ChemEquil::nAtoms(), Phase::speciesName(), and Phase::temperature().

Referenced by ChemEquil::equilibrate(), ChemEquil::setInitialMoles(), and ChemEquil::setToEquilState().

double calcEmoles	(	thermo_t &	s,
		vector_fp &	x,
		const double &	n_t,
		const vector_fp &	Xmol_i_calc,
		vector_fp &	eMolesCalc,
		vector_fp &	n_i_calc,
		double	pressureConst
	)

protected

Given a vector of dimensionless element abundances, this routine calculates the moles of the elements and the moles of the species.

Input

x[m] = current dimensionless element potentials..

Definition at line 1250 of file ChemEquil.cpp.

References DATA_PTR, ThermoPhase::getActivityCoefficients(), ChemEquil::m_muSS_RT, ChemEquil::nAtoms(), Phase::setMoleFractions(), and ThermoPhase::setPressure().

Referenced by ChemEquil::estimateEP_Brinkley().

Member Data Documentation

EquilOpt options

Options controlling how the calculation is carried out.

See Also: EquilOptions

Definition at line 127 of file ChemEquil.h.

Referenced by Cantera::equilibrate(), ChemEquil::equilibrate(), ChemEquil::estimateEP_Brinkley(), and Cantera::vcs_equilibrate().

thermo_t* m_phase

protected

Pointer to the ThermoPhase object used to initialize this object.

This ThermoPhase object must be compatible with the ThermoPhase objects input from the equilibrate function. Currently, this means that the 2 ThermoPhases have to have consist of the same species and elements.

Definition at line 140 of file ChemEquil.h.

Referenced by ChemEquil::initialize().

size_t m_nComponents

protected

This is equal to the rank of the stoichiometric coefficient matrix when it is computed.

It's initialized to m_mm.

Definition at line 188 of file ChemEquil.h.

Referenced by ChemEquil::equilibrate(), ChemEquil::equilResidual(), ChemEquil::estimateElementPotentials(), ChemEquil::estimateEP_Brinkley(), ChemEquil::initialize(), and ChemEquil::setInitialMoles().

vector_fp m_molefractions

protected

Current value of the mole fractions in the single phase.

-> length = m_kk.

Definition at line 196 of file ChemEquil.h.

Referenced by ChemEquil::initialize(), and ChemEquil::update().

vector_fp m_lambda

protected

Current value of the dimensional element potentials -> length = m_mm.

Definition at line 201 of file ChemEquil.h.

Referenced by ChemEquil::equilibrate(), and ChemEquil::initialize().

size_t m_eloc

protected

Index of the element id corresponding to the electric charge of each species.

Equal to -1 if there is no such element id.

Definition at line 227 of file ChemEquil.h.

Referenced by ChemEquil::equilibrate(), ChemEquil::equilResidual(), ChemEquil::estimateEP_Brinkley(), and ChemEquil::initialize().

vector_fp m_muSS_RT

protected

Dimensionless values of the gibbs free energy for the standard state of each species, at the temperature and pressure of the solution (the star standard state).

Definition at line 238 of file ChemEquil.h.

Referenced by ChemEquil::calcEmoles(), ChemEquil::estimateEP_Brinkley(), and ChemEquil::initialize().

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

Public Member Functions

Public Attributes

Protected Member Functions

Protected Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Input

Member Data Documentation