1558 lines
53 KiB
C++
1558 lines
53 KiB
C++
/**
|
|
* @file vcs_MultiPhaseEquil.cpp
|
|
* Driver routine for the VCSnonideal equilibrium solver package
|
|
*/
|
|
/*
|
|
* Copyright (2006) Sandia Corporation. Under the terms of
|
|
* Contract DE-AC04-94AL85000 with Sandia Corporation, the
|
|
* U.S. Government retains certain rights in this software.
|
|
*/
|
|
|
|
#include "cantera/equil/vcs_MultiPhaseEquil.h"
|
|
#include "cantera/equil/vcs_prob.h"
|
|
#include "cantera/equil/vcs_internal.h"
|
|
#include "cantera/equil/vcs_VolPhase.h"
|
|
#include "cantera/equil/vcs_species_thermo.h"
|
|
|
|
#include "cantera/equil/vcs_solve.h"
|
|
|
|
#include "cantera/thermo/mix_defs.h"
|
|
#include "cantera/base/clockWC.h"
|
|
#include "cantera/base/stringUtils.h"
|
|
#include "cantera/thermo/speciesThermoTypes.h"
|
|
#include "cantera/thermo/IdealSolidSolnPhase.h"
|
|
#include "cantera/thermo/IdealMolalSoln.h"
|
|
#include "cantera/equil/ChemEquil.h"
|
|
#include "cantera/equil/vcs_SpeciesProperties.h"
|
|
|
|
#include <cstdio>
|
|
|
|
using namespace Cantera;
|
|
using namespace std;
|
|
|
|
namespace VCSnonideal
|
|
{
|
|
vcs_MultiPhaseEquil::vcs_MultiPhaseEquil() :
|
|
m_vprob(0),
|
|
m_mix(0),
|
|
m_printLvl(0),
|
|
m_vsolvePtr(0)
|
|
{
|
|
}
|
|
|
|
vcs_MultiPhaseEquil::vcs_MultiPhaseEquil(Cantera::MultiPhase* mix, int printLvl) :
|
|
m_vprob(0),
|
|
m_mix(0),
|
|
m_printLvl(printLvl),
|
|
m_vsolvePtr(0)
|
|
{
|
|
/*
|
|
* Create a VCS_PROB object that describes the equilibrium problem.
|
|
* The constructor just mallocs the necessary objects and sizes them.
|
|
*/
|
|
m_vprob = new VCS_PROB(mix->nSpecies(),
|
|
mix->nElements(),
|
|
mix->nPhases());
|
|
m_mix = mix;
|
|
m_vprob->m_printLvl = m_printLvl;
|
|
/*
|
|
* Work out the details of the VCS_VPROB construction and
|
|
* Transfer the current problem to VCS_PROB object
|
|
*/
|
|
int res = vcs_Cantera_to_vprob(mix, m_vprob);
|
|
if (res != 0) {
|
|
plogf("problems\n");
|
|
}
|
|
}
|
|
|
|
vcs_MultiPhaseEquil::~vcs_MultiPhaseEquil()
|
|
{
|
|
delete m_vprob;
|
|
m_vprob = 0;
|
|
delete m_vsolvePtr;
|
|
m_vsolvePtr = 0;
|
|
}
|
|
|
|
int vcs_MultiPhaseEquil::equilibrate_TV(int XY, doublereal xtarget,
|
|
int estimateEquil,
|
|
int printLvl, doublereal err,
|
|
int maxsteps, int loglevel)
|
|
{
|
|
// doublereal dt = 1.0e3;
|
|
doublereal Vtarget = m_mix->volume();
|
|
doublereal dVdP;
|
|
if ((XY != TV) && (XY != HV) && (XY != UV) && (XY != SV)) {
|
|
throw CanteraError("vcs_MultiPhaseEquil::equilibrate_TV",
|
|
"Wrong XY flag:" + int2str(XY));
|
|
}
|
|
int maxiter = 100;
|
|
int iSuccess = 0;
|
|
int innerXY;
|
|
double Pnow;
|
|
if (XY == TV) {
|
|
m_mix->setTemperature(xtarget);
|
|
}
|
|
double Pnew;
|
|
int strt = estimateEquil;
|
|
double P1 = 0.0;
|
|
double V1 = 0.0;
|
|
double V2 = 0.0;
|
|
double P2 = 0.0;
|
|
doublereal Tlow = 0.5 * m_mix->minTemp();
|
|
doublereal Thigh = 2.0 * m_mix->maxTemp();
|
|
doublereal Vnow, Verr;
|
|
int printLvlSub = std::max(0, printLvl - 1);
|
|
for (int n = 0; n < maxiter; n++) {
|
|
Pnow = m_mix->pressure();
|
|
|
|
switch (XY) {
|
|
case TV:
|
|
iSuccess = equilibrate_TP(strt, printLvlSub, err, maxsteps, loglevel);
|
|
break;
|
|
case HV:
|
|
innerXY = HP;
|
|
iSuccess = equilibrate_HP(xtarget, innerXY, Tlow, Thigh, strt,
|
|
printLvlSub, err, maxsteps, loglevel);
|
|
break;
|
|
case UV:
|
|
innerXY = UP;
|
|
iSuccess = equilibrate_HP(xtarget, innerXY, Tlow, Thigh, strt,
|
|
printLvlSub, err, maxsteps, loglevel);
|
|
break;
|
|
case SV:
|
|
innerXY = SP;
|
|
iSuccess = equilibrate_SP(xtarget, Tlow, Thigh, strt,
|
|
printLvlSub, err, maxsteps, loglevel);
|
|
break;
|
|
default:
|
|
break;
|
|
}
|
|
strt = false;
|
|
Vnow = m_mix->volume();
|
|
if (n == 0) {
|
|
V2 = Vnow;
|
|
P2 = Pnow;
|
|
} else if (n == 1) {
|
|
V1 = Vnow;
|
|
P1 = Pnow;
|
|
} else {
|
|
P2 = P1;
|
|
V2 = V1;
|
|
P1 = Pnow;
|
|
V1 = Vnow;
|
|
}
|
|
|
|
Verr = fabs((Vtarget - Vnow)/Vtarget);
|
|
|
|
if (Verr < err) {
|
|
goto done;
|
|
}
|
|
// find dV/dP
|
|
if (n > 1) {
|
|
dVdP = (V2 - V1) / (P2 - P1);
|
|
if (dVdP == 0.0) {
|
|
throw CanteraError("vcs_MultiPhase::equilibrate_TV",
|
|
"dVdP == 0.0");
|
|
} else {
|
|
Pnew = Pnow + (Vtarget - Vnow) / dVdP;
|
|
if (Pnew < 0.2 * Pnow) {
|
|
Pnew = 0.2 * Pnow;
|
|
}
|
|
if (Pnew > 3.0 * Pnow) {
|
|
Pnew = 3.0 * Pnow;
|
|
}
|
|
}
|
|
|
|
} else {
|
|
m_mix->setPressure(Pnow*1.01);
|
|
dVdP = (m_mix->volume() - Vnow)/(0.01*Pnow);
|
|
Pnew = Pnow + 0.5*(Vtarget - Vnow)/dVdP;
|
|
if (Pnew < 0.5* Pnow) {
|
|
Pnew = 0.5 * Pnow;
|
|
}
|
|
if (Pnew > 1.7 * Pnow) {
|
|
Pnew = 1.7 * Pnow;
|
|
}
|
|
|
|
}
|
|
m_mix->setPressure(Pnew);
|
|
}
|
|
throw CanteraError("vcs_MultiPhase::equilibrate_TV",
|
|
"No convergence for V");
|
|
|
|
done:
|
|
;
|
|
return iSuccess;
|
|
}
|
|
|
|
int vcs_MultiPhaseEquil::equilibrate_HP(doublereal Htarget,
|
|
int XY, double Tlow, double Thigh,
|
|
int estimateEquil,
|
|
int printLvl, doublereal err,
|
|
int maxsteps, int loglevel)
|
|
{
|
|
int maxiter = 100;
|
|
int iSuccess;
|
|
if (XY != HP && XY != UP) {
|
|
throw CanteraError("vcs_MultiPhaseEquil::equilibrate_HP",
|
|
"Wrong XP" + XY);
|
|
}
|
|
int strt = estimateEquil;
|
|
|
|
// Lower bound on T. This will change as we progress in the calculation
|
|
if (Tlow <= 0.0) {
|
|
Tlow = 0.5 * m_mix->minTemp();
|
|
}
|
|
// Upper bound on T. This will change as we progress in the calculation
|
|
if (Thigh <= 0.0 || Thigh > 1.0E6) {
|
|
Thigh = 2.0 * m_mix->maxTemp();
|
|
}
|
|
|
|
doublereal cpb = 1.0, dT, dTa, dTmax, Tnew;
|
|
doublereal Hnow;
|
|
doublereal Hlow = Undef;
|
|
doublereal Hhigh = Undef;
|
|
doublereal Herr, HConvErr;
|
|
doublereal Tnow = m_mix->temperature();
|
|
int printLvlSub = std::max(printLvl - 1, 0);
|
|
|
|
for (int n = 0; n < maxiter; n++) {
|
|
|
|
|
|
// start with a loose error tolerance, but tighten it as we get
|
|
// close to the final temperature
|
|
try {
|
|
Tnow = m_mix->temperature();
|
|
iSuccess = equilibrate_TP(strt, printLvlSub, err, maxsteps, loglevel);
|
|
strt = 0;
|
|
if (XY == UP) {
|
|
Hnow = m_mix->IntEnergy();
|
|
} else {
|
|
Hnow = m_mix->enthalpy();
|
|
}
|
|
double pmoles[10];
|
|
pmoles[0] = m_mix->phaseMoles(0);
|
|
double Tmoles = pmoles[0];
|
|
double HperMole = Hnow/Tmoles;
|
|
if (printLvl > 0) {
|
|
plogf("T = %g, Hnow = %g ,Tmoles = %g, HperMole = %g",
|
|
Tnow, Hnow, Tmoles, HperMole);
|
|
plogendl();
|
|
}
|
|
|
|
// the equilibrium enthalpy monotonically increases with T;
|
|
// if the current value is below the target, then we know the
|
|
// current temperature is too low. Set the lower bounds.
|
|
|
|
|
|
if (Hnow < Htarget) {
|
|
if (Tnow > Tlow) {
|
|
Tlow = Tnow;
|
|
Hlow = Hnow;
|
|
}
|
|
}
|
|
// the current enthalpy is greater than the target; therefore the
|
|
// current temperature is too high. Set the high bounds.
|
|
else {
|
|
if (Tnow < Thigh) {
|
|
Thigh = Tnow;
|
|
Hhigh = Hnow;
|
|
}
|
|
}
|
|
if (Hlow != Undef && Hhigh != Undef) {
|
|
cpb = (Hhigh - Hlow)/(Thigh - Tlow);
|
|
dT = (Htarget - Hnow)/cpb;
|
|
dTa = fabs(dT);
|
|
dTmax = 0.5*fabs(Thigh - Tlow);
|
|
if (dTa > dTmax) {
|
|
dT *= dTmax/dTa;
|
|
}
|
|
} else {
|
|
Tnew = sqrt(Tlow*Thigh);
|
|
dT = Tnew - Tnow;
|
|
if (dT < -200.) {
|
|
dT = -200;
|
|
}
|
|
if (dT > 200.) {
|
|
dT = 200.;
|
|
}
|
|
}
|
|
double acpb = std::max(fabs(cpb), 1.0E-6);
|
|
double denom = std::max(fabs(Htarget), acpb);
|
|
Herr = Htarget - Hnow;
|
|
HConvErr = fabs((Herr)/denom);
|
|
if (printLvl > 0) {
|
|
plogf(" equilibrate_HP: It = %d, Tcurr = %g Hcurr = %g, Htarget = %g\n",
|
|
n, Tnow, Hnow, Htarget);
|
|
plogf(" H rel error = %g, cp = %g, HConvErr = %g\n",
|
|
Herr, cpb, HConvErr);
|
|
}
|
|
|
|
if (HConvErr < err) { // || dTa < 1.0e-4) {
|
|
if (printLvl > 0) {
|
|
plogf(" equilibrate_HP: CONVERGENCE: Hfinal = %g Tfinal = %g, Its = %d \n",
|
|
Hnow, Tnow, n);
|
|
plogf(" H rel error = %g, cp = %g, HConvErr = %g\n",
|
|
Herr, cpb, HConvErr);
|
|
}
|
|
goto done;
|
|
}
|
|
Tnew = Tnow + dT;
|
|
if (Tnew < 0.0) {
|
|
Tnew = 0.5*Tnow;
|
|
}
|
|
m_mix->setTemperature(Tnew);
|
|
|
|
} catch (CanteraError err) {
|
|
if (!estimateEquil) {
|
|
strt = -1;
|
|
} else {
|
|
Tnew = 0.5*(Tnow + Thigh);
|
|
if (fabs(Tnew - Tnow) < 1.0) {
|
|
Tnew = Tnow + 1.0;
|
|
}
|
|
m_mix->setTemperature(Tnew);
|
|
}
|
|
}
|
|
|
|
}
|
|
throw CanteraError("MultiPhase::equilibrate_HP",
|
|
"No convergence for T");
|
|
done:
|
|
;
|
|
return iSuccess;
|
|
}
|
|
|
|
int vcs_MultiPhaseEquil::equilibrate_SP(doublereal Starget,
|
|
double Tlow, double Thigh,
|
|
int estimateEquil,
|
|
int printLvl, doublereal err,
|
|
int maxsteps, int loglevel)
|
|
{
|
|
int maxiter = 100;
|
|
int strt = estimateEquil;
|
|
|
|
// Lower bound on T. This will change as we progress in the calculation
|
|
if (Tlow <= 0.0) {
|
|
Tlow = 0.5 * m_mix->minTemp();
|
|
}
|
|
// Upper bound on T. This will change as we progress in the calculation
|
|
if (Thigh <= 0.0 || Thigh > 1.0E6) {
|
|
Thigh = 2.0 * m_mix->maxTemp();
|
|
}
|
|
|
|
doublereal cpb = 1.0, dT, dTa, dTmax, Tnew;
|
|
doublereal Snow;
|
|
doublereal Slow = Undef;
|
|
doublereal Shigh = Undef;
|
|
doublereal Serr, SConvErr;
|
|
doublereal Tnow = m_mix->temperature();
|
|
if (Tnow < Tlow) {
|
|
Tlow = Tnow;
|
|
}
|
|
if (Tnow > Thigh) {
|
|
Thigh = Tnow;
|
|
}
|
|
int printLvlSub = std::max(printLvl - 1, 0);
|
|
|
|
for (int n = 0; n < maxiter; n++) {
|
|
|
|
// start with a loose error tolerance, but tighten it as we get
|
|
// close to the final temperature
|
|
try {
|
|
Tnow = m_mix->temperature();
|
|
int iSuccess = equilibrate_TP(strt, printLvlSub, err, maxsteps, loglevel);
|
|
strt = 0;
|
|
Snow = m_mix->entropy();
|
|
double pmoles[10];
|
|
pmoles[0] = m_mix->phaseMoles(0);
|
|
double Tmoles = pmoles[0];
|
|
double SperMole = Snow/Tmoles;
|
|
plogf("T = %g, Snow = %g ,Tmoles = %g, SperMole = %g\n",
|
|
Tnow, Snow, Tmoles, SperMole);
|
|
|
|
// the equilibrium entropy monotonically increases with T;
|
|
// if the current value is below the target, then we know the
|
|
// current temperature is too low. Set the lower bounds to the
|
|
// current condition.
|
|
if (Snow < Starget) {
|
|
if (Tnow > Tlow) {
|
|
Tlow = Tnow;
|
|
Slow = Snow;
|
|
} else {
|
|
if (Slow > Starget) {
|
|
if (Snow < Slow) {
|
|
Thigh = Tlow;
|
|
Shigh = Slow;
|
|
Tlow = Tnow;
|
|
Slow = Snow;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
// the current enthalpy is greater than the target; therefore the
|
|
// current temperature is too high. Set the high bounds.
|
|
else {
|
|
if (Tnow < Thigh) {
|
|
Thigh = Tnow;
|
|
Shigh = Snow;
|
|
}
|
|
}
|
|
if (Slow != Undef && Shigh != Undef) {
|
|
cpb = (Shigh - Slow)/(Thigh - Tlow);
|
|
dT = (Starget - Snow)/cpb;
|
|
Tnew = Tnow + dT;
|
|
dTa = fabs(dT);
|
|
dTmax = 0.5*fabs(Thigh - Tlow);
|
|
if (Tnew > Thigh || Tnew < Tlow) {
|
|
dTmax = 1.5*fabs(Thigh - Tlow);
|
|
}
|
|
dTmax = std::min(dTmax, 300.);
|
|
if (dTa > dTmax) {
|
|
dT *= dTmax/dTa;
|
|
}
|
|
} else {
|
|
Tnew = sqrt(Tlow*Thigh);
|
|
dT = Tnew - Tnow;
|
|
}
|
|
|
|
double acpb = std::max(fabs(cpb), 1.0E-6);
|
|
double denom = std::max(fabs(Starget), acpb);
|
|
Serr = Starget - Snow;
|
|
SConvErr = fabs((Serr)/denom);
|
|
if (printLvl > 0) {
|
|
plogf(" equilibrate_SP: It = %d, Tcurr = %g Scurr = %g, Starget = %g\n",
|
|
n, Tnow, Snow, Starget);
|
|
plogf(" S rel error = %g, cp = %g, SConvErr = %g\n",
|
|
Serr, cpb, SConvErr);
|
|
}
|
|
|
|
if (SConvErr < err) { // || dTa < 1.0e-4) {
|
|
if (printLvl > 0) {
|
|
plogf(" equilibrate_SP: CONVERGENCE: Sfinal = %g Tfinal = %g, Its = %d \n",
|
|
Snow, Tnow, n);
|
|
plogf(" S rel error = %g, cp = %g, HConvErr = %g\n",
|
|
Serr, cpb, SConvErr);
|
|
}
|
|
return iSuccess;
|
|
}
|
|
Tnew = Tnow + dT;
|
|
if (Tnew < 0.0) {
|
|
Tnew = 0.5*Tnow;
|
|
}
|
|
m_mix->setTemperature(Tnew);
|
|
|
|
} catch (CanteraError err) {
|
|
if (!estimateEquil) {
|
|
strt = -1;
|
|
} else {
|
|
Tnew = 0.5*(Tnow + Thigh);
|
|
if (fabs(Tnew - Tnow) < 1.0) {
|
|
Tnew = Tnow + 1.0;
|
|
}
|
|
m_mix->setTemperature(Tnew);
|
|
}
|
|
}
|
|
|
|
}
|
|
throw CanteraError("MultiPhase::equilibrate_SP",
|
|
"No convergence for T");
|
|
}
|
|
|
|
int vcs_MultiPhaseEquil::equilibrate(int XY, int estimateEquil,
|
|
int printLvl, doublereal err,
|
|
int maxsteps, int loglevel)
|
|
{
|
|
int iSuccess;
|
|
doublereal xtarget;
|
|
if (XY == TP) {
|
|
iSuccess = equilibrate_TP(estimateEquil, printLvl, err, maxsteps, loglevel);
|
|
} else if (XY == HP || XY == UP) {
|
|
if (XY == HP) {
|
|
xtarget = m_mix->enthalpy();
|
|
} else {
|
|
xtarget = m_mix->IntEnergy();
|
|
}
|
|
double Tlow = 0.5 * m_mix->minTemp();
|
|
double Thigh = 2.0 * m_mix->maxTemp();
|
|
iSuccess = equilibrate_HP(xtarget, XY, Tlow, Thigh,
|
|
estimateEquil, printLvl, err, maxsteps, loglevel);
|
|
} else if (XY == SP) {
|
|
xtarget = m_mix->entropy();
|
|
double Tlow = 0.5 * m_mix->minTemp();
|
|
double Thigh = 2.0 * m_mix->maxTemp();
|
|
iSuccess = equilibrate_SP(xtarget, Tlow, Thigh,
|
|
estimateEquil, printLvl, err, maxsteps, loglevel);
|
|
|
|
} else if (XY == TV) {
|
|
xtarget = m_mix->temperature();
|
|
iSuccess = equilibrate_TV(XY, xtarget,
|
|
estimateEquil, printLvl, err, maxsteps, loglevel);
|
|
} else if (XY == HV) {
|
|
xtarget = m_mix->enthalpy();
|
|
iSuccess = equilibrate_TV(XY, xtarget,
|
|
estimateEquil, printLvl, err, maxsteps, loglevel);
|
|
} else if (XY == UV) {
|
|
xtarget = m_mix->IntEnergy();
|
|
iSuccess = equilibrate_TV(XY, xtarget,
|
|
estimateEquil, printLvl, err, maxsteps, loglevel);
|
|
} else if (XY == SV) {
|
|
xtarget = m_mix->entropy();
|
|
iSuccess = equilibrate_TV(XY, xtarget, estimateEquil,
|
|
printLvl, err, maxsteps, loglevel);
|
|
} else {
|
|
throw CanteraError(" vcs_MultiPhaseEquil::equilibrate",
|
|
"Unsupported Option");
|
|
}
|
|
return iSuccess;
|
|
}
|
|
|
|
int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil,
|
|
int printLvl, doublereal err,
|
|
int maxsteps, int loglevel)
|
|
{
|
|
int maxit = maxsteps;
|
|
clockWC tickTock;
|
|
|
|
if (m_vprob == 0) {
|
|
m_vprob = new VCS_PROB(m_mix->nSpecies(),
|
|
m_mix->nElements(),
|
|
m_mix->nPhases());
|
|
}
|
|
m_printLvl = printLvl;
|
|
m_vprob->m_printLvl = printLvl;
|
|
|
|
|
|
/*
|
|
* Extract the current state information
|
|
* from the MultiPhase object and
|
|
* Transfer it to VCS_PROB object.
|
|
*/
|
|
int res = vcs_Cantera_update_vprob(m_mix, m_vprob);
|
|
if (res != 0) {
|
|
plogf("problems\n");
|
|
}
|
|
|
|
|
|
// Set the estimation technique
|
|
if (estimateEquil) {
|
|
m_vprob->iest = estimateEquil;
|
|
} else {
|
|
m_vprob->iest = 0;
|
|
}
|
|
|
|
// Check obvious bounds on the temperature and pressure
|
|
// NOTE, we may want to do more here with the real bounds
|
|
// given by the ThermoPhase objects.
|
|
double T = m_mix->temperature();
|
|
if (T <= 0.0) {
|
|
throw CanteraError("vcs_MultiPhaseEquil::equilibrate",
|
|
"Temperature less than zero on input");
|
|
}
|
|
double pres = m_mix->pressure();
|
|
if (pres <= 0.0) {
|
|
throw CanteraError("vcs_MultiPhaseEquil::equilibrate",
|
|
"Pressure less than zero on input");
|
|
}
|
|
|
|
/*
|
|
* Print out the problem specification from the point of
|
|
* view of the vprob object.
|
|
*/
|
|
m_vprob->prob_report(m_printLvl);
|
|
|
|
/*
|
|
* Call the thermo Program
|
|
*/
|
|
int ip1 = m_printLvl;
|
|
int ipr = std::max(0, m_printLvl-1);
|
|
if (m_printLvl >= 3) {
|
|
ip1 = m_printLvl - 2;
|
|
} else {
|
|
ip1 = 0;
|
|
}
|
|
if (!m_vsolvePtr) {
|
|
m_vsolvePtr = new VCS_SOLVE();
|
|
}
|
|
int iSuccess = m_vsolvePtr->vcs(m_vprob, 0, ipr, ip1, maxit);
|
|
|
|
/*
|
|
* Transfer the information back to the MultiPhase object.
|
|
* Note we don't just call setMoles, because some multispecies
|
|
* solution phases may be zeroed out, and that would cause a problem
|
|
* for that routine. Also, the mole fractions of such zeroed out
|
|
* phases actually contain information about likely reemergent
|
|
* states.
|
|
*/
|
|
m_mix->uploadMoleFractionsFromPhases();
|
|
size_t kGlob = 0;
|
|
for (size_t ip = 0; ip < m_vprob->NPhase; ip++) {
|
|
double phaseMole = 0.0;
|
|
Cantera::ThermoPhase& tref = m_mix->phase(ip);
|
|
for (size_t k = 0; k < tref.nSpecies(); k++, kGlob++) {
|
|
phaseMole += m_vprob->w[kGlob];
|
|
}
|
|
//phaseMole *= 1.0E-3;
|
|
m_mix->setPhaseMoles(ip, phaseMole);
|
|
}
|
|
|
|
double te = tickTock.secondsWC();
|
|
if (printLvl > 0) {
|
|
plogf("\n Results from vcs:\n");
|
|
if (iSuccess != 0) {
|
|
plogf("\nVCS FAILED TO CONVERGE!\n");
|
|
}
|
|
plogf("\n");
|
|
plogf("Temperature = %g Kelvin\n", m_vprob->T);
|
|
plogf("Pressure = %g Pa\n", m_vprob->PresPA);
|
|
plogf("\n");
|
|
plogf("----------------------------------------"
|
|
"---------------------\n");
|
|
plogf(" Name Mole_Number");
|
|
if (m_vprob->m_VCS_UnitsFormat == VCS_UNITS_MKS) {
|
|
plogf("(kmol)");
|
|
} else {
|
|
plogf("(gmol)");
|
|
}
|
|
plogf(" Mole_Fraction Chem_Potential");
|
|
if (m_vprob->m_VCS_UnitsFormat == VCS_UNITS_KCALMOL) {
|
|
plogf(" (kcal/mol)\n");
|
|
} else if (m_vprob->m_VCS_UnitsFormat == VCS_UNITS_UNITLESS) {
|
|
plogf(" (Dimensionless)\n");
|
|
} else if (m_vprob->m_VCS_UnitsFormat == VCS_UNITS_KJMOL) {
|
|
plogf(" (kJ/mol)\n");
|
|
} else if (m_vprob->m_VCS_UnitsFormat == VCS_UNITS_KELVIN) {
|
|
plogf(" (Kelvin)\n");
|
|
} else if (m_vprob->m_VCS_UnitsFormat == VCS_UNITS_MKS) {
|
|
plogf(" (J/kmol)\n");
|
|
}
|
|
plogf("--------------------------------------------------"
|
|
"-----------\n");
|
|
for (size_t i = 0; i < m_vprob->nspecies; i++) {
|
|
plogf("%-12s", m_vprob->SpName[i].c_str());
|
|
if (m_vprob->SpeciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
|
|
plogf(" %15.3e %15.3e ", 0.0, m_vprob->mf[i]);
|
|
plogf("%15.3e\n", m_vprob->m_gibbsSpecies[i]);
|
|
} else {
|
|
plogf(" %15.3e %15.3e ", m_vprob->w[i], m_vprob->mf[i]);
|
|
if (m_vprob->w[i] <= 0.0) {
|
|
size_t iph = m_vprob->PhaseID[i];
|
|
vcs_VolPhase* VPhase = m_vprob->VPhaseList[iph];
|
|
if (VPhase->nSpecies() > 1) {
|
|
plogf(" -1.000e+300\n");
|
|
} else {
|
|
plogf("%15.3e\n", m_vprob->m_gibbsSpecies[i]);
|
|
}
|
|
} else {
|
|
plogf("%15.3e\n", m_vprob->m_gibbsSpecies[i]);
|
|
}
|
|
}
|
|
}
|
|
plogf("------------------------------------------"
|
|
"-------------------\n");
|
|
if (printLvl > 2) {
|
|
if (m_vsolvePtr->m_timing_print_lvl > 0) {
|
|
plogf("Total time = %12.6e seconds\n", te);
|
|
}
|
|
}
|
|
}
|
|
return iSuccess;
|
|
}
|
|
|
|
void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile)
|
|
{
|
|
size_t k;
|
|
size_t istart;
|
|
size_t nSpecies;
|
|
|
|
double vol = 0.0;
|
|
string sName;
|
|
size_t nphase = m_vprob->NPhase;
|
|
|
|
FILE* FP = fopen(reportFile.c_str(), "w");
|
|
if (!FP) {
|
|
plogf("Failure to open file\n");
|
|
exit(EXIT_FAILURE);
|
|
}
|
|
double Temp = m_mix->temperature();
|
|
double pres = m_mix->pressure();
|
|
double* mf = VCS_DATA_PTR(m_vprob->mf);
|
|
#ifdef DEBUG_MODE
|
|
double* fe = VCS_DATA_PTR(m_vprob->m_gibbsSpecies);
|
|
#endif
|
|
std::vector<double> VolPM;
|
|
std::vector<double> activity;
|
|
std::vector<double> ac;
|
|
std::vector<double> mu;
|
|
std::vector<double> mu0;
|
|
std::vector<double> molalities;
|
|
|
|
|
|
vol = 0.0;
|
|
for (size_t iphase = 0; iphase < nphase; iphase++) {
|
|
istart = m_mix->speciesIndex(0, iphase);
|
|
Cantera::ThermoPhase& tref = m_mix->phase(iphase);
|
|
nSpecies = tref.nSpecies();
|
|
VolPM.resize(nSpecies, 0.0);
|
|
tref.getPartialMolarVolumes(VCS_DATA_PTR(VolPM));
|
|
vcs_VolPhase* volP = m_vprob->VPhaseList[iphase];
|
|
|
|
double TMolesPhase = volP->totalMoles();
|
|
double VolPhaseVolumes = 0.0;
|
|
for (k = 0; k < nSpecies; k++) {
|
|
VolPhaseVolumes += VolPM[k] * mf[istart + k];
|
|
}
|
|
VolPhaseVolumes *= TMolesPhase;
|
|
vol += VolPhaseVolumes;
|
|
}
|
|
|
|
fprintf(FP,"--------------------- VCS_MULTIPHASE_EQUIL FINAL REPORT"
|
|
" -----------------------------\n");
|
|
fprintf(FP,"Temperature = %11.5g kelvin\n", Temp);
|
|
fprintf(FP,"Pressure = %11.5g Pascal\n", pres);
|
|
fprintf(FP,"Total Volume = %11.5g m**3\n", vol);
|
|
fprintf(FP,"Number Basis optimizations = %d\n", m_vprob->m_NumBasisOptimizations);
|
|
fprintf(FP,"Number VCS iterations = %d\n", m_vprob->m_Iterations);
|
|
|
|
for (size_t iphase = 0; iphase < nphase; iphase++) {
|
|
istart = m_mix->speciesIndex(0, iphase);
|
|
Cantera::ThermoPhase& tref = m_mix->phase(iphase);
|
|
Cantera::ThermoPhase* tp = &tref;
|
|
string phaseName = tref.name();
|
|
vcs_VolPhase* volP = m_vprob->VPhaseList[iphase];
|
|
double TMolesPhase = volP->totalMoles();
|
|
//AssertTrace(TMolesPhase == m_mix->phaseMoles(iphase));
|
|
nSpecies = tref.nSpecies();
|
|
activity.resize(nSpecies, 0.0);
|
|
ac.resize(nSpecies, 0.0);
|
|
|
|
mu0.resize(nSpecies, 0.0);
|
|
mu.resize(nSpecies, 0.0);
|
|
VolPM.resize(nSpecies, 0.0);
|
|
molalities.resize(nSpecies, 0.0);
|
|
|
|
int actConvention = tp->activityConvention();
|
|
tp->getActivities(VCS_DATA_PTR(activity));
|
|
tp->getActivityCoefficients(VCS_DATA_PTR(ac));
|
|
tp->getStandardChemPotentials(VCS_DATA_PTR(mu0));
|
|
|
|
tp->getPartialMolarVolumes(VCS_DATA_PTR(VolPM));
|
|
tp->getChemPotentials(VCS_DATA_PTR(mu));
|
|
double VolPhaseVolumes = 0.0;
|
|
for (k = 0; k < nSpecies; k++) {
|
|
VolPhaseVolumes += VolPM[k] * mf[istart + k];
|
|
}
|
|
VolPhaseVolumes *= TMolesPhase;
|
|
vol += VolPhaseVolumes;
|
|
|
|
|
|
if (actConvention == 1) {
|
|
MolalityVPSSTP* mTP = static_cast<MolalityVPSSTP*>(tp);
|
|
mTP->getMolalities(VCS_DATA_PTR(molalities));
|
|
tp->getChemPotentials(VCS_DATA_PTR(mu));
|
|
|
|
if (iphase == 0) {
|
|
fprintf(FP," Name, Phase, PhaseMoles, Mole_Fract, "
|
|
"Molalities, ActCoeff, Activity,"
|
|
"ChemPot_SS0, ChemPot, mole_num, PMVol, Phase_Volume\n");
|
|
|
|
fprintf(FP," , , (kmol), , "
|
|
" , , ,"
|
|
" (J/kmol), (J/kmol), (kmol), (m**3/kmol), (m**3)\n");
|
|
}
|
|
for (k = 0; k < nSpecies; k++) {
|
|
sName = tp->speciesName(k);
|
|
fprintf(FP,"%12s, %11s, %11.3e, %11.3e, %11.3e, %11.3e, %11.3e,"
|
|
"%11.3e, %11.3e, %11.3e, %11.3e, %11.3e\n",
|
|
sName.c_str(),
|
|
phaseName.c_str(), TMolesPhase,
|
|
mf[istart + k], molalities[k], ac[k], activity[k],
|
|
mu0[k]*1.0E-6, mu[k]*1.0E-6,
|
|
mf[istart + k] * TMolesPhase,
|
|
VolPM[k], VolPhaseVolumes);
|
|
}
|
|
|
|
} else {
|
|
if (iphase == 0) {
|
|
fprintf(FP," Name, Phase, PhaseMoles, Mole_Fract, "
|
|
"Molalities, ActCoeff, Activity,"
|
|
" ChemPotSS0, ChemPot, mole_num, PMVol, Phase_Volume\n");
|
|
|
|
fprintf(FP," , , (kmol), , "
|
|
" , , ,"
|
|
" (J/kmol), (J/kmol), (kmol), (m**3/kmol), (m**3)\n");
|
|
}
|
|
for (k = 0; k < nSpecies; k++) {
|
|
molalities[k] = 0.0;
|
|
}
|
|
for (k = 0; k < nSpecies; k++) {
|
|
sName = tp->speciesName(k);
|
|
fprintf(FP,"%12s, %11s, %11.3e, %11.3e, %11.3e, %11.3e, %11.3e, "
|
|
"%11.3e, %11.3e,% 11.3e, %11.3e, %11.3e\n",
|
|
sName.c_str(),
|
|
phaseName.c_str(), TMolesPhase,
|
|
mf[istart + k], molalities[k], ac[k],
|
|
activity[k], mu0[k]*1.0E-6, mu[k]*1.0E-6,
|
|
mf[istart + k] * TMolesPhase,
|
|
VolPM[k], VolPhaseVolumes);
|
|
}
|
|
}
|
|
|
|
#ifdef DEBUG_MODE
|
|
/*
|
|
* Check consistency: These should be equal
|
|
*/
|
|
tp->getChemPotentials(fe+istart);
|
|
for (k = 0; k < nSpecies; k++) {
|
|
if (!vcs_doubleEqual(fe[istart+k], mu[k])) {
|
|
fprintf(FP,"ERROR: incompatibility!\n");
|
|
fclose(FP);
|
|
plogf("ERROR: incompatibility!\n");
|
|
exit(EXIT_FAILURE);
|
|
}
|
|
}
|
|
#endif
|
|
|
|
}
|
|
fclose(FP);
|
|
}
|
|
|
|
//! print char repeatedly to log file
|
|
/*!
|
|
* @param letter letter to be repeated
|
|
* @param num Number of times repeated
|
|
*/
|
|
static void print_char(const char letter, const int num)
|
|
{
|
|
for (int i = 0; i < num; i++) {
|
|
plogf("%c", letter);
|
|
}
|
|
}
|
|
|
|
/*
|
|
* HKM -> Work on transferring the current value of the voltages into the
|
|
* equilibrium problem.
|
|
*/
|
|
int vcs_Cantera_to_vprob(Cantera::MultiPhase* mphase,
|
|
VCSnonideal::VCS_PROB* vprob)
|
|
{
|
|
VCS_SPECIES_THERMO* ts_ptr = 0;
|
|
|
|
/*
|
|
* Calculate the total number of species and phases in the problem
|
|
*/
|
|
size_t totNumPhases = mphase->nPhases();
|
|
size_t totNumSpecies = mphase->nSpecies();
|
|
|
|
// Problem type has yet to be worked out.
|
|
vprob->prob_type = 0;
|
|
vprob->nspecies = totNumSpecies;
|
|
vprob->ne = 0;
|
|
vprob->NPhase = totNumPhases;
|
|
vprob->m_VCS_UnitsFormat = VCS_UNITS_MKS;
|
|
// Set the initial estimate to a machine generated estimate for now
|
|
// We will work out the details later.
|
|
vprob->iest = -1;
|
|
vprob->T = mphase->temperature();
|
|
vprob->PresPA = mphase->pressure();
|
|
vprob->Vol = mphase->volume();
|
|
vprob->Title = "MultiPhase Object";
|
|
|
|
Cantera::ThermoPhase* tPhase = 0;
|
|
|
|
bool gasPhase;
|
|
int printLvl = vprob->m_printLvl;
|
|
|
|
/*
|
|
* Loop over the phases, transferring pertinent information
|
|
*/
|
|
int kT = 0;
|
|
for (size_t iphase = 0; iphase < totNumPhases; iphase++) {
|
|
|
|
/*
|
|
* Get the thermophase object - assume volume phase
|
|
*/
|
|
tPhase = &(mphase->phase(iphase));
|
|
size_t nelem = tPhase->nElements();
|
|
|
|
/*
|
|
* Query Cantera for the equation of state type of the
|
|
* current phase.
|
|
*/
|
|
int eos = tPhase->eosType();
|
|
if (eos == cIdealGas) {
|
|
gasPhase = true;
|
|
} else {
|
|
gasPhase = false;
|
|
}
|
|
|
|
/*
|
|
* Find out the number of species in the phase
|
|
*/
|
|
size_t nSpPhase = tPhase->nSpecies();
|
|
/*
|
|
* Find out the name of the phase
|
|
*/
|
|
string phaseName = tPhase->name();
|
|
|
|
/*
|
|
* Call the basic vcs_VolPhase creation routine.
|
|
* Properties set here:
|
|
* ->PhaseNum = phase number in the thermo problem
|
|
* ->GasPhase = Boolean indicating whether it is a gas phase
|
|
* ->NumSpecies = number of species in the phase
|
|
* ->TMolesInert = Inerts in the phase = 0.0 for cantera
|
|
* ->PhaseName = Name of the phase
|
|
*/
|
|
|
|
vcs_VolPhase* VolPhase = vprob->VPhaseList[iphase];
|
|
VolPhase->resize(iphase, nSpPhase, nelem, phaseName.c_str(), 0.0);
|
|
VolPhase->m_gasPhase = gasPhase;
|
|
/*
|
|
* Tell the vcs_VolPhase pointer about cantera
|
|
*/
|
|
VolPhase->p_VCS_UnitsFormat = vprob->m_VCS_UnitsFormat;
|
|
VolPhase->setPtrThermoPhase(tPhase);
|
|
VolPhase->setTotalMoles(0.0);
|
|
/*
|
|
* Set the electric potential of the volume phase from the
|
|
* ThermoPhase object's value.
|
|
*/
|
|
VolPhase->setElectricPotential(tPhase->electricPotential());
|
|
/*
|
|
* Query the ThermoPhase object to find out what convention
|
|
* it uses for the specification of activity and Standard State.
|
|
*/
|
|
VolPhase->p_activityConvention = tPhase->activityConvention();
|
|
/*
|
|
* Assign the value of eqn of state
|
|
* -> Handle conflicts here.
|
|
*/
|
|
switch (eos) {
|
|
case cIdealGas:
|
|
VolPhase->m_eqnState = VCS_EOS_IDEAL_GAS;
|
|
break;
|
|
case cIncompressible:
|
|
VolPhase->m_eqnState = VCS_EOS_CONSTANT;
|
|
break;
|
|
case cSurf:
|
|
throw CanteraError("VCSnonideal", "cSurf not handled yet.");
|
|
case cStoichSubstance:
|
|
VolPhase->m_eqnState = VCS_EOS_STOICH_SUB;
|
|
break;
|
|
case cPureFluid:
|
|
if (printLvl > 1) {
|
|
plogf("cPureFluid not recognized yet by VCSnonideal\n");
|
|
}
|
|
break;
|
|
case cEdge:
|
|
throw CanteraError("VCSnonideal", "cEdge not handled yet.");
|
|
case cIdealSolidSolnPhase0:
|
|
case cIdealSolidSolnPhase1:
|
|
case cIdealSolidSolnPhase2:
|
|
VolPhase->m_eqnState = VCS_EOS_IDEAL_SOLN;
|
|
break;
|
|
default:
|
|
if (printLvl > 1) {
|
|
plogf("Unknown Cantera EOS to VCSnonideal: %d\n", eos);
|
|
}
|
|
VolPhase->m_eqnState = VCS_EOS_UNK_CANTERA;
|
|
if (!VolPhase->usingCanteraCalls()) {
|
|
plogf("vcs functions asked for, but unimplemented\n");
|
|
exit(EXIT_FAILURE);
|
|
}
|
|
break;
|
|
}
|
|
|
|
/*
|
|
* Transfer all of the element information from the
|
|
* ThermoPhase object to the vcs_VolPhase object.
|
|
* Also decide whether we need a new charge neutrality
|
|
* element in the phase to enforce a charge neutrality
|
|
* constraint.
|
|
* We also decide whether this is a single species phase
|
|
* with the voltage being the independent variable setting
|
|
* the chemical potential of the electrons.
|
|
*/
|
|
VolPhase->transferElementsFM(tPhase);
|
|
|
|
/*
|
|
* Combine the element information in the vcs_VolPhase
|
|
* object into the vprob object.
|
|
*/
|
|
vprob->addPhaseElements(VolPhase);
|
|
|
|
VolPhase->setState_TP(vprob->T, vprob->PresPA);
|
|
vector<double> muPhase(tPhase->nSpecies(),0.0);
|
|
tPhase->getChemPotentials(&muPhase[0]);
|
|
double tMoles = 0.0;
|
|
/*
|
|
* Loop through each species in the current phase
|
|
*/
|
|
for (size_t k = 0; k < nSpPhase; k++) {
|
|
/*
|
|
* Obtain the molecular weight of the species from the
|
|
* ThermoPhase object
|
|
*/
|
|
vprob->WtSpecies[kT] = tPhase->molecularWeight(k);
|
|
|
|
/*
|
|
* Obtain the charges of the species from the
|
|
* ThermoPhase object
|
|
*/
|
|
vprob->Charge[kT] = tPhase->charge(k);
|
|
|
|
/*
|
|
* Set the phaseid of the species
|
|
*/
|
|
vprob->PhaseID[kT] = iphase;
|
|
|
|
/*
|
|
* Transfer the Species name
|
|
*/
|
|
string stmp = mphase->speciesName(kT);
|
|
vprob->SpName[kT] = stmp;
|
|
/*
|
|
* Transfer the type of unknown
|
|
*/
|
|
vprob->SpeciesUnknownType[kT] = VolPhase->speciesUnknownType(k);
|
|
|
|
if (vprob->SpeciesUnknownType[kT] == VCS_SPECIES_TYPE_MOLNUM) {
|
|
/*
|
|
* Set the initial number of kmoles of the species
|
|
* and the mole fraction vector
|
|
*/
|
|
vprob->w[kT] = mphase->speciesMoles(kT);
|
|
tMoles += vprob->w[kT];
|
|
vprob->mf[kT] = mphase->moleFraction(kT);
|
|
} else if (vprob->SpeciesUnknownType[kT] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
|
|
vprob->w[kT] = tPhase->electricPotential();
|
|
vprob->mf[kT] = mphase->moleFraction(kT);
|
|
} else {
|
|
throw CanteraError(" vcs_Cantera_to_vprob() ERROR",
|
|
"Unknown species type: " + int2str(vprob->SpeciesUnknownType[kT]));
|
|
}
|
|
|
|
|
|
/*
|
|
* transfer chemical potential vector
|
|
*/
|
|
vprob->m_gibbsSpecies[kT] = muPhase[k];
|
|
/*
|
|
* Transfer the species information from the
|
|
* volPhase structure to the VPROB structure
|
|
* This includes:
|
|
* FormulaMatrix[][]
|
|
* VolPhase->IndSpecies[]
|
|
*/
|
|
vprob->addOnePhaseSpecies(VolPhase, k, kT);
|
|
|
|
/*
|
|
* Get a pointer to the thermo object
|
|
*/
|
|
ts_ptr = vprob->SpeciesThermo[kT];
|
|
/*
|
|
* Fill in the vcs_SpeciesProperty structure
|
|
*/
|
|
vcs_SpeciesProperties* sProp = VolPhase->speciesProperty(k);
|
|
sProp->NumElements = vprob->ne;
|
|
sProp->SpName = vprob->SpName[kT];
|
|
sProp->SpeciesThermo = ts_ptr;
|
|
sProp->WtSpecies = tPhase->molecularWeight(k);
|
|
sProp->FormulaMatrixCol.resize(vprob->ne, 0.0);
|
|
for (size_t e = 0; e < vprob->ne; e++) {
|
|
sProp->FormulaMatrixCol[e] = vprob->FormulaMatrix[e][kT];
|
|
}
|
|
sProp->Charge = tPhase->charge(k);
|
|
sProp->SurfaceSpecies = false;
|
|
sProp->VolPM = 0.0;
|
|
|
|
/*
|
|
* Transfer the thermo specification of the species
|
|
* vprob->SpeciesThermo[]
|
|
*/
|
|
ts_ptr->UseCanteraCalls = VolPhase->usingCanteraCalls();
|
|
ts_ptr->m_VCS_UnitsFormat = VolPhase->p_VCS_UnitsFormat;
|
|
/*
|
|
* Add lookback connectivity into the thermo object first
|
|
*/
|
|
ts_ptr->IndexPhase = iphase;
|
|
ts_ptr->IndexSpeciesPhase = k;
|
|
ts_ptr->OwningPhase = VolPhase;
|
|
/*
|
|
* get a reference to the Cantera species thermo.
|
|
*/
|
|
SpeciesThermo& sp = tPhase->speciesThermo();
|
|
|
|
int spType = sp.reportType(k);
|
|
|
|
if (spType == SIMPLE) {
|
|
double c[4];
|
|
double minTemp, maxTemp, refPressure;
|
|
sp.reportParams(k, spType, c, minTemp, maxTemp, refPressure);
|
|
ts_ptr->SS0_Model = VCS_SS0_CONSTANT;
|
|
ts_ptr->SS0_T0 = c[0];
|
|
ts_ptr->SS0_H0 = c[1];
|
|
ts_ptr->SS0_S0 = c[2];
|
|
ts_ptr->SS0_Cp0 = c[3];
|
|
if (gasPhase) {
|
|
ts_ptr->SSStar_Model = VCS_SSSTAR_IDEAL_GAS;
|
|
ts_ptr->SSStar_Vol_Model = VCS_SSVOL_IDEALGAS;
|
|
} else {
|
|
ts_ptr->SSStar_Model = VCS_SSSTAR_CONSTANT;
|
|
ts_ptr->SSStar_Vol_Model = VCS_SSVOL_CONSTANT;
|
|
}
|
|
ts_ptr->Activity_Coeff_Model = VCS_AC_CONSTANT;
|
|
ts_ptr->Activity_Coeff_Params = NULL;
|
|
} else {
|
|
if (vprob->m_printLvl > 2) {
|
|
plogf("vcs_Cantera_convert: Species Type %d not known \n",
|
|
spType);
|
|
}
|
|
ts_ptr->SS0_Model = VCS_SS0_NOTHANDLED;
|
|
ts_ptr->SSStar_Model = VCS_SSSTAR_NOTHANDLED;
|
|
if (!(ts_ptr->UseCanteraCalls)) {
|
|
plogf("Cantera calls not being used -> exiting\n");
|
|
exit(EXIT_FAILURE);
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Transfer the Volume Information -> NEEDS WORK
|
|
*/
|
|
if (gasPhase) {
|
|
ts_ptr->SSStar_Vol_Model = VCS_SSVOL_IDEALGAS;
|
|
ts_ptr->SSStar_Vol_Params = NULL;
|
|
ts_ptr->SSStar_Vol0 = 82.05 * 273.15 / 1.0;
|
|
|
|
} else {
|
|
std::vector<double> phaseTermCoeff(nSpPhase, 0.0);
|
|
int nCoeff;
|
|
tPhase->getParameters(nCoeff, VCS_DATA_PTR(phaseTermCoeff));
|
|
ts_ptr->SSStar_Vol_Model = VCS_SSVOL_CONSTANT;
|
|
ts_ptr->SSStar_Vol0 = phaseTermCoeff[k];
|
|
}
|
|
kT++;
|
|
}
|
|
|
|
/*
|
|
* Now go back through the species in the phase and assign
|
|
* a valid mole fraction to all phases, even if the initial
|
|
* estimate of the total number of moles is zero.
|
|
*/
|
|
if (tMoles > 0.0) {
|
|
for (size_t k = 0; k < nSpPhase; k++) {
|
|
size_t kTa = VolPhase->spGlobalIndexVCS(k);
|
|
vprob->mf[kTa] = vprob->w[kTa] / tMoles;
|
|
}
|
|
} else {
|
|
/*
|
|
* Perhaps, we could do a more sophisticated treatment below.
|
|
* But, will start with this.
|
|
*/
|
|
for (size_t k = 0; k < nSpPhase; k++) {
|
|
size_t kTa = VolPhase->spGlobalIndexVCS(k);
|
|
vprob->mf[kTa]= 1.0 / (double) nSpPhase;
|
|
}
|
|
}
|
|
|
|
VolPhase->setMolesFromVCS(VCS_STATECALC_OLD, VCS_DATA_PTR(vprob->w));
|
|
/*
|
|
* Now, calculate a sample naught gibbs free energy calculation
|
|
* at the specified temperature.
|
|
*/
|
|
double R = vcsUtil_gasConstant(vprob->m_VCS_UnitsFormat);
|
|
for (size_t k = 0; k < nSpPhase; k++) {
|
|
vcs_SpeciesProperties* sProp = VolPhase->speciesProperty(k);
|
|
ts_ptr = sProp->SpeciesThermo;
|
|
ts_ptr->SS0_feSave = VolPhase->G0_calc_one(k)/ R;
|
|
ts_ptr->SS0_TSave = vprob->T;
|
|
}
|
|
|
|
}
|
|
|
|
/*
|
|
* Transfer initial element abundances to the vprob object.
|
|
* We have to find the mapping index from one to the other
|
|
*
|
|
*/
|
|
vprob->gai.resize(vprob->ne, 0.0);
|
|
vprob->set_gai();
|
|
|
|
/*
|
|
* Printout the species information: PhaseID's and mole nums
|
|
*/
|
|
if (vprob->m_printLvl > 1) {
|
|
plogf("\n");
|
|
print_char('=', 80);
|
|
plogf("\n");
|
|
print_char('=', 16);
|
|
plogf(" Cantera_to_vprob: START OF PROBLEM STATEMENT ");
|
|
print_char('=', 20);
|
|
plogf("\n");
|
|
print_char('=', 80);
|
|
plogf("\n");
|
|
plogf(" Phase IDs of species\n");
|
|
plogf(" species phaseID phaseName ");
|
|
plogf(" Initial_Estimated_kMols\n");
|
|
for (size_t i = 0; i < vprob->nspecies; i++) {
|
|
size_t iphase = vprob->PhaseID[i];
|
|
|
|
vcs_VolPhase* VolPhase = vprob->VPhaseList[iphase];
|
|
plogf("%16s %5d %16s", vprob->SpName[i].c_str(), iphase,
|
|
VolPhase->PhaseName.c_str());
|
|
if (vprob->SpeciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
|
|
plogf(" Volts = %-10.5g\n", vprob->w[i]);
|
|
} else {
|
|
plogf(" %-10.5g\n", vprob->w[i]);
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Printout of the Phase structure information
|
|
*/
|
|
plogf("\n");
|
|
print_char('-', 80);
|
|
plogf("\n");
|
|
plogf(" Information about phases\n");
|
|
plogf(" PhaseName PhaseNum SingSpec GasPhase EqnState NumSpec");
|
|
plogf(" TMolesInert Tmoles(kmol)\n");
|
|
|
|
for (size_t iphase = 0; iphase < vprob->NPhase; iphase++) {
|
|
vcs_VolPhase* VolPhase = vprob->VPhaseList[iphase];
|
|
std::string sEOS = string16_EOSType(VolPhase->m_eqnState);
|
|
plogf("%16s %5d %5d %8d %16s %8d %16e ", VolPhase->PhaseName.c_str(),
|
|
VolPhase->VP_ID_, VolPhase->m_singleSpecies,
|
|
VolPhase->m_gasPhase, sEOS.c_str(),
|
|
VolPhase->nSpecies(), VolPhase->totalMolesInert());
|
|
plogf("%16e\n", VolPhase->totalMoles());
|
|
}
|
|
|
|
plogf("\n");
|
|
print_char('=', 80);
|
|
plogf("\n");
|
|
print_char('=', 16);
|
|
plogf(" Cantera_to_vprob: END OF PROBLEM STATEMENT ");
|
|
print_char('=', 20);
|
|
plogf("\n");
|
|
print_char('=', 80);
|
|
plogf("\n\n");
|
|
}
|
|
|
|
return VCS_SUCCESS;
|
|
}
|
|
|
|
int vcs_Cantera_update_vprob(Cantera::MultiPhase* mphase,
|
|
VCSnonideal::VCS_PROB* vprob)
|
|
{
|
|
size_t totNumPhases = mphase->nPhases();
|
|
size_t kT = 0;
|
|
std::vector<double> tmpMoles;
|
|
// Problem type has yet to be worked out.
|
|
vprob->prob_type = 0;
|
|
// Whether we have an estimate or not gets overwritten on
|
|
// the call to the equilibrium solver.
|
|
vprob->iest = -1;
|
|
vprob->T = mphase->temperature();
|
|
vprob->PresPA = mphase->pressure();
|
|
vprob->Vol = mphase->volume();
|
|
Cantera::ThermoPhase* tPhase = 0;
|
|
|
|
for (size_t iphase = 0; iphase < totNumPhases; iphase++) {
|
|
tPhase = &(mphase->phase(iphase));
|
|
vcs_VolPhase* volPhase = vprob->VPhaseList[iphase];
|
|
/*
|
|
* Set the electric potential of the volume phase from the
|
|
* ThermoPhase object's value.
|
|
*/
|
|
volPhase->setElectricPotential(tPhase->electricPotential());
|
|
|
|
volPhase->setState_TP(vprob->T, vprob->PresPA);
|
|
vector<double> muPhase(tPhase->nSpecies(),0.0);
|
|
tPhase->getChemPotentials(&muPhase[0]);
|
|
/*
|
|
* Loop through each species in the current phase
|
|
*/
|
|
size_t nSpPhase = tPhase->nSpecies();
|
|
// volPhase->TMoles = 0.0;
|
|
tmpMoles.resize(nSpPhase);
|
|
for (size_t k = 0; k < nSpPhase; k++) {
|
|
tmpMoles[k] = mphase->speciesMoles(kT);
|
|
vprob->w[kT] = mphase->speciesMoles(kT);
|
|
vprob->mf[kT] = mphase->moleFraction(kT);
|
|
|
|
/*
|
|
* transfer chemical potential vector
|
|
*/
|
|
vprob->m_gibbsSpecies[kT] = muPhase[k];
|
|
|
|
kT++;
|
|
}
|
|
if (volPhase->phiVarIndex() != npos) {
|
|
size_t kphi = volPhase->phiVarIndex();
|
|
size_t kglob = volPhase->spGlobalIndexVCS(kphi);
|
|
vprob->w[kglob] = tPhase->electricPotential();
|
|
}
|
|
volPhase->setMolesFromVCS(VCS_STATECALC_OLD, VCS_DATA_PTR(vprob->w));
|
|
if ((nSpPhase == 1) && (volPhase->phiVarIndex() == 0)) {
|
|
volPhase->setExistence(VCS_PHASE_EXIST_ALWAYS);
|
|
} else if (volPhase->totalMoles() > 0.0) {
|
|
volPhase->setExistence(VCS_PHASE_EXIST_YES);
|
|
} else {
|
|
volPhase->setExistence(VCS_PHASE_EXIST_NO);
|
|
}
|
|
|
|
}
|
|
/*
|
|
* Transfer initial element abundances to the vprob object.
|
|
* Put them in the front of the object. There may be
|
|
* more constraints than there are elements. But, we
|
|
* know the element abundances are in the front of the
|
|
* vector.
|
|
*/
|
|
vprob->set_gai();
|
|
|
|
/*
|
|
* Printout the species information: PhaseID's and mole nums
|
|
*/
|
|
if (vprob->m_printLvl > 1) {
|
|
plogf("\n");
|
|
print_char('=', 80);
|
|
plogf("\n");
|
|
print_char('=', 20);
|
|
plogf(" Cantera_to_vprob: START OF PROBLEM STATEMENT ");
|
|
print_char('=', 20);
|
|
plogf("\n");
|
|
print_char('=', 80);
|
|
plogf("\n\n");
|
|
plogf(" Phase IDs of species\n");
|
|
plogf(" species phaseID phaseName ");
|
|
plogf(" Initial_Estimated_kMols\n");
|
|
for (size_t i = 0; i < vprob->nspecies; i++) {
|
|
size_t iphase = vprob->PhaseID[i];
|
|
|
|
vcs_VolPhase* VolPhase = vprob->VPhaseList[iphase];
|
|
plogf("%16s %5d %16s", vprob->SpName[i].c_str(), iphase,
|
|
VolPhase->PhaseName.c_str());
|
|
if (vprob->SpeciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
|
|
plogf(" Volts = %-10.5g\n", vprob->w[i]);
|
|
} else {
|
|
plogf(" %-10.5g\n", vprob->w[i]);
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Printout of the Phase structure information
|
|
*/
|
|
plogf("\n");
|
|
print_char('-', 80);
|
|
plogf("\n");
|
|
plogf(" Information about phases\n");
|
|
plogf(" PhaseName PhaseNum SingSpec GasPhase EqnState NumSpec");
|
|
plogf(" TMolesInert Tmoles(kmol)\n");
|
|
|
|
for (size_t iphase = 0; iphase < vprob->NPhase; iphase++) {
|
|
vcs_VolPhase* VolPhase = vprob->VPhaseList[iphase];
|
|
std::string sEOS = string16_EOSType(VolPhase->m_eqnState);
|
|
plogf("%16s %5d %5d %8d %16s %8d %16e ", VolPhase->PhaseName.c_str(),
|
|
VolPhase->VP_ID_, VolPhase->m_singleSpecies,
|
|
VolPhase->m_gasPhase, sEOS.c_str(),
|
|
VolPhase->nSpecies(), VolPhase->totalMolesInert());
|
|
plogf("%16e\n", VolPhase->totalMoles());
|
|
}
|
|
|
|
plogf("\n");
|
|
print_char('=', 80);
|
|
plogf("\n");
|
|
print_char('=', 20);
|
|
plogf(" Cantera_to_vprob: END OF PROBLEM STATEMENT ");
|
|
print_char('=', 20);
|
|
plogf("\n");
|
|
print_char('=', 80);
|
|
plogf("\n\n");
|
|
}
|
|
|
|
return VCS_SUCCESS;
|
|
}
|
|
|
|
void vcs_MultiPhaseEquil::getStoichVector(size_t rxn, Cantera::vector_fp& nu)
|
|
{
|
|
size_t nsp = m_vsolvePtr->m_numSpeciesTot;
|
|
nu.resize(nsp, 0.0);
|
|
for (size_t i = 0; i < nsp; i++) {
|
|
nu[i] = 0.0;
|
|
}
|
|
size_t nc = numComponents();
|
|
// scMatrix [nrxn][ncomp]
|
|
const DoubleStarStar& scMatrix = m_vsolvePtr->m_stoichCoeffRxnMatrix;
|
|
const std::vector<size_t>& indSpecies = m_vsolvePtr->m_speciesMapIndex;
|
|
if (rxn > nsp - nc) {
|
|
return;
|
|
}
|
|
size_t j = indSpecies[rxn + nc];
|
|
nu[j] = 1.0;
|
|
for (size_t kc = 0; kc < nc; kc++) {
|
|
j = indSpecies[kc];
|
|
nu[j] = scMatrix[rxn][kc];
|
|
}
|
|
|
|
}
|
|
|
|
size_t vcs_MultiPhaseEquil::numComponents() const
|
|
{
|
|
size_t nc = npos;
|
|
if (m_vsolvePtr) {
|
|
nc = m_vsolvePtr->m_numComponents;
|
|
}
|
|
return nc;
|
|
}
|
|
|
|
size_t vcs_MultiPhaseEquil::numElemConstraints() const
|
|
{
|
|
size_t nec = npos;
|
|
if (m_vsolvePtr) {
|
|
nec = m_vsolvePtr->m_numElemConstraints;
|
|
}
|
|
return nec;
|
|
}
|
|
|
|
size_t vcs_MultiPhaseEquil::component(size_t m) const
|
|
{
|
|
size_t nc = numComponents();
|
|
if (m < nc) {
|
|
return m_vsolvePtr->m_speciesMapIndex[m];
|
|
} else {
|
|
return npos;
|
|
}
|
|
}
|
|
|
|
int vcs_MultiPhaseEquil::determine_PhaseStability(int iph, double& funcStab, int printLvl, int loglevel)
|
|
{
|
|
clockWC tickTock;
|
|
size_t nsp = m_mix->nSpecies();
|
|
size_t nel = m_mix->nElements();
|
|
size_t nph = m_mix->nPhases();
|
|
if (m_vprob == 0) {
|
|
m_vprob = new VCS_PROB(nsp, nel, nph);
|
|
}
|
|
m_printLvl = printLvl;
|
|
m_vprob->m_printLvl = printLvl;
|
|
|
|
/*
|
|
* Extract the current state information
|
|
* from the MultiPhase object and
|
|
* Transfer it to VCS_PROB object.
|
|
*/
|
|
int res = vcs_Cantera_update_vprob(m_mix, m_vprob);
|
|
if (res != 0) {
|
|
plogf("problems\n");
|
|
}
|
|
|
|
|
|
|
|
// Check obvious bounds on the temperature and pressure
|
|
// NOTE, we may want to do more here with the real bounds
|
|
// given by the ThermoPhase objects.
|
|
double T = m_mix->temperature();
|
|
if (T <= 0.0) {
|
|
throw CanteraError("vcs_MultiPhaseEquil::determine_PhaseStability",
|
|
"Temperature less than zero on input");
|
|
}
|
|
double pres = m_mix->pressure();
|
|
if (pres <= 0.0) {
|
|
throw CanteraError("vcs_MultiPhaseEquil::determine_PhaseStability",
|
|
"Pressure less than zero on input");
|
|
}
|
|
|
|
/*
|
|
* Print out the problem specification from the point of
|
|
* view of the vprob object.
|
|
*/
|
|
m_vprob->prob_report(m_printLvl);
|
|
|
|
/*
|
|
* Call the thermo Program
|
|
*/
|
|
if (!m_vsolvePtr) {
|
|
m_vsolvePtr = new VCS_SOLVE();
|
|
}
|
|
int iStable = m_vsolvePtr->vcs_PS(m_vprob, iph, printLvl, funcStab);
|
|
|
|
/*
|
|
* Transfer the information back to the MultiPhase object.
|
|
* Note we don't just call setMoles, because some multispecies
|
|
* solution phases may be zeroed out, and that would cause a problem
|
|
* for that routine. Also, the mole fractions of such zeroed out
|
|
* phases actually contain information about likely reemergent
|
|
* states.
|
|
*/
|
|
m_mix->uploadMoleFractionsFromPhases();
|
|
// for (int i = 0; i < m_vprob->nspecies; i++) {
|
|
// plogf("%d %15.3e\n", m_vprob->m_gibbsSpecies[i]);
|
|
//}
|
|
m_mix->getChemPotentials(DATA_PTR(m_vprob->m_gibbsSpecies));
|
|
//for (int i = 0; i < m_vprob->nspecies; i++) {
|
|
// plogf("%d %15.3e\n", m_vprob->m_gibbsSpecies[i]);
|
|
//}
|
|
|
|
double te = tickTock.secondsWC();
|
|
if (printLvl > 0) {
|
|
plogf("\n Results from vcs_PS:\n");
|
|
|
|
plogf("\n");
|
|
plogf("Temperature = %g Kelvin\n", m_vprob->T);
|
|
plogf("Pressure = %g Pa\n", m_vprob->PresPA);
|
|
std::string sss = m_mix->phaseName(iph);
|
|
if (iStable) {
|
|
plogf("Phase %d named %s is stable, function value = %g > 0\n", iph, sss.c_str(), funcStab);
|
|
} else {
|
|
plogf("Phase %d named %s is not stable + function value = %g < 0\n", iph, sss.c_str(), funcStab);
|
|
}
|
|
plogf("\n");
|
|
plogf("----------------------------------------"
|
|
"---------------------\n");
|
|
plogf(" Name Mole_Number");
|
|
if (m_vprob->m_VCS_UnitsFormat == VCS_UNITS_MKS) {
|
|
plogf("(kmol)");
|
|
} else {
|
|
plogf("(gmol)");
|
|
}
|
|
plogf(" Mole_Fraction Chem_Potential");
|
|
if (m_vprob->m_VCS_UnitsFormat == VCS_UNITS_KCALMOL) {
|
|
plogf(" (kcal/mol)\n");
|
|
} else if (m_vprob->m_VCS_UnitsFormat == VCS_UNITS_UNITLESS) {
|
|
plogf(" (Dimensionless)\n");
|
|
} else if (m_vprob->m_VCS_UnitsFormat == VCS_UNITS_KJMOL) {
|
|
plogf(" (kJ/mol)\n");
|
|
} else if (m_vprob->m_VCS_UnitsFormat == VCS_UNITS_KELVIN) {
|
|
plogf(" (Kelvin)\n");
|
|
} else if (m_vprob->m_VCS_UnitsFormat == VCS_UNITS_MKS) {
|
|
plogf(" (J/kmol)\n");
|
|
}
|
|
plogf("-------------------------------------------------------------\n");
|
|
for (size_t i = 0; i < m_vprob->nspecies; i++) {
|
|
plogf("%-12s", m_vprob->SpName[i].c_str());
|
|
if (m_vprob->SpeciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
|
|
plogf(" %15.3e %15.3e ", 0.0, m_vprob->mf[i]);
|
|
plogf("%15.3e\n", m_vprob->m_gibbsSpecies[i]);
|
|
} else {
|
|
plogf(" %15.3e %15.3e ", m_vprob->w[i], m_vprob->mf[i]);
|
|
if (m_vprob->w[i] <= 0.0) {
|
|
plogf("%15.3e\n", m_vprob->m_gibbsSpecies[i]);
|
|
} else {
|
|
plogf("%15.3e\n", m_vprob->m_gibbsSpecies[i]);
|
|
}
|
|
}
|
|
}
|
|
plogf("------------------------------------------"
|
|
"-------------------\n");
|
|
if (printLvl > 2) {
|
|
if (m_vsolvePtr->m_timing_print_lvl > 0) {
|
|
plogf("Total time = %12.6e seconds\n", te);
|
|
}
|
|
}
|
|
}
|
|
return iStable;
|
|
}
|
|
|
|
}
|