diff --git a/Cantera/src/equil/vcs_Gibbs.cpp b/Cantera/src/equil/vcs_Gibbs.cpp index 3b0b00b52..5fa7ea891 100644 --- a/Cantera/src/equil/vcs_Gibbs.cpp +++ b/Cantera/src/equil/vcs_Gibbs.cpp @@ -45,11 +45,7 @@ namespace VCSnonideal { g += TPhInertMoles[iph] * log(TPhInertMoles[iph] / tPhMoles[iph]); if (Vphase->GasPhase) { - if (m_VCS_UnitsFormat == 3) { - g += TPhInertMoles[iph] * log(m_pressure/(Cantera::OneAtm)); - } else { - g += TPhInertMoles[iph] * log(m_pressure); - } + g += TPhInertMoles[iph] * log(m_pressurePA/(1.01325E5)); } } } @@ -86,11 +82,7 @@ namespace VCSnonideal { g += TPhInertMoles[iphase] * log(TPhInertMoles[iphase] / phaseMols); vcs_VolPhase *Vphase = VPhaseList[iphase]; if (Vphase->GasPhase == iphase) { - if (m_VCS_UnitsFormat == 3) { - g += TPhInertMoles[iphase] * log(m_pressure/(Cantera::OneAtm)); - } else { - g += TPhInertMoles[iphase] * log(m_pressure); - } + g += TPhInertMoles[iphase] * log(m_pressurePA/1.01325E5); } } diff --git a/Cantera/src/equil/vcs_MultiPhaseEquil.cpp b/Cantera/src/equil/vcs_MultiPhaseEquil.cpp index c1bbd7967..7a4cedad4 100644 --- a/Cantera/src/equil/vcs_MultiPhaseEquil.cpp +++ b/Cantera/src/equil/vcs_MultiPhaseEquil.cpp @@ -668,12 +668,7 @@ namespace VCSnonideal { } plogf("\n"); plogf("Temperature = %g Kelvin\n", m_vprob->T); - plogf("Pressure = %g ", m_vprob->Pres); - if (m_vprob->m_VCS_UnitsFormat == VCS_UNITS_MKS) { - plogf("Pa\n"); - } else { - plogf("atm\n"); - } + plogf("Pressure = %g Pa\n", m_vprob->PresPA); plogf("\n"); plogf("----------------------------------------" "---------------------\n"); @@ -1099,7 +1094,7 @@ namespace VCSnonideal { // We will work out the details later. vprob->iest = -1; vprob->T = mphase->temperature(); - vprob->Pres = mphase->pressure(); + vprob->PresPA = mphase->pressure(); vprob->Vol = mphase->volume(); vprob->Title = "MultiPhase Object"; @@ -1225,7 +1220,7 @@ namespace VCSnonideal { */ vprob->addPhaseElements(VolPhase); - VolPhase->setState_TP(vprob->T, vprob->Pres); + VolPhase->setState_TP(vprob->T, vprob->PresPA); vector muPhase(tPhase->nSpecies(),0.0); tPhase->getChemPotentials(&muPhase[0]); double tMoles = 0.0; @@ -1426,14 +1421,14 @@ namespace VCSnonideal { print_char('=', 80); plogf("\n"); plogf(" Phase IDs of species\n"); plogf(" species phaseID phaseName "); - plogf(" Initial_Estimated_gMols\n"); + plogf(" Initial_Estimated_kMols\n"); for (int i = 0; i < vprob->nspecies; i++) { int iphase = vprob->PhaseID[i]; vcs_VolPhase *VolPhase = vprob->VPhaseList[iphase]; plogf("%16s %5d %16s", vprob->SpName[i].c_str(), iphase, VolPhase->PhaseName.c_str()); - plogf(" %-10.5g\n", 1.0E3 * vprob->w[i]); + plogf(" %-10.5g\n", vprob->w[i]); } /* @@ -1442,7 +1437,7 @@ namespace VCSnonideal { plogf("\n"); print_char('-', 80); plogf("\n"); plogf(" Information about phases\n"); plogf(" PhaseName PhaseNum SingSpec GasPhase EqnState NumSpec"); - plogf(" TMolesInert Tmoles(gmol)\n"); + plogf(" TMolesInert Tmoles(kmol)\n"); for (int iphase = 0; iphase < vprob->NPhase; iphase++) { vcs_VolPhase *VolPhase = vprob->VPhaseList[iphase]; @@ -1450,8 +1445,8 @@ namespace VCSnonideal { plogf("%16s %5d %5d %8d %16s %8d %16e ", VolPhase->PhaseName.c_str(), VolPhase->VP_ID, VolPhase->SingleSpecies, VolPhase->GasPhase, sEOS.c_str(), - VolPhase->NVolSpecies, VolPhase->TMolesInert * 1.0E3); - plogf("%16e\n", VolPhase->TotalMoles() * 1.0E3); + VolPhase->NVolSpecies, VolPhase->TMolesInert ); + plogf("%16e\n", VolPhase->TotalMoles()); } plogf("\n"); print_char('=', 80); plogf("\n"); @@ -1479,7 +1474,7 @@ namespace VCSnonideal { // the call to the equilibrium solver. vprob->iest = -1; vprob->T = mphase->temperature(); - vprob->Pres = mphase->pressure(); + vprob->PresPA = mphase->pressure(); vprob->Vol = mphase->volume(); Cantera::ThermoPhase *tPhase = 0; @@ -1492,7 +1487,7 @@ namespace VCSnonideal { */ volPhase->setElectricPotential(tPhase->electricPotential()); - volPhase->setState_TP(vprob->T, vprob->Pres); + volPhase->setState_TP(vprob->T, vprob->PresPA); vector muPhase(tPhase->nSpecies(),0.0); tPhase->getChemPotentials(&muPhase[0]); /* @@ -1544,14 +1539,14 @@ namespace VCSnonideal { print_char('=', 80); plogf("\n\n"); plogf(" Phase IDs of species\n"); plogf(" species phaseID phaseName "); - plogf(" Initial_Estimated_gMols\n"); + plogf(" Initial_Estimated_kMols\n"); for (int i = 0; i < vprob->nspecies; i++) { int iphase = vprob->PhaseID[i]; vcs_VolPhase *VolPhase = vprob->VPhaseList[iphase]; plogf("%16s %5d %16s", vprob->SpName[i].c_str(), iphase, VolPhase->PhaseName.c_str()); - plogf(" %-10.5g\n", 1.0E3 * vprob->w[i]); + plogf(" %-10.5g\n", vprob->w[i]); } /* @@ -1560,7 +1555,7 @@ namespace VCSnonideal { plogf("\n"); print_char('-', 80); plogf("\n"); plogf(" Information about phases\n"); plogf(" PhaseName PhaseNum SingSpec GasPhase EqnState NumSpec"); - plogf(" TMolesInert Tmoles(gmol)\n"); + plogf(" TMolesInert Tmoles(kmol)\n"); for (int iphase = 0; iphase < vprob->NPhase; iphase++) { vcs_VolPhase *VolPhase = vprob->VPhaseList[iphase]; @@ -1568,8 +1563,8 @@ namespace VCSnonideal { plogf("%16s %5d %5d %8d %16s %8d %16e ", VolPhase->PhaseName.c_str(), VolPhase->VP_ID, VolPhase->SingleSpecies, VolPhase->GasPhase, sEOS.c_str(), - VolPhase->NVolSpecies, VolPhase->TMolesInert * 1.0E3); - plogf("%16e\n", VolPhase->TotalMoles() * 1.0E3); + VolPhase->NVolSpecies, VolPhase->TMolesInert ); + plogf("%16e\n", VolPhase->TotalMoles() ); } plogf("\n"); print_char('=', 80); plogf("\n"); diff --git a/Cantera/src/equil/vcs_TP.cpp b/Cantera/src/equil/vcs_TP.cpp index 772a4f726..19c3df783 100644 --- a/Cantera/src/equil/vcs_TP.cpp +++ b/Cantera/src/equil/vcs_TP.cpp @@ -36,7 +36,7 @@ namespace VCSnonideal { * ip1 = 1 -> Print intermediate results. * maxit -> Maximum number of iterations for the algorithm * T = Temperature (Kelvin) - * pres = Pressure (units given by if__ variable) + * pres = Pressure (pascal) * * Return Codes * ------------------ @@ -55,7 +55,7 @@ namespace VCSnonideal { * Store the temperature and pressure in the private global variables */ m_temperature = T_arg; - m_pressure = pres_arg; + m_pressurePA = pres_arg; /* * Evaluate the standard state free energies * at the current temperatures and pressures. @@ -113,7 +113,7 @@ namespace VCSnonideal { * 0 -> don't report on anything * IP1 = 1 -> Print intermediate results. * T = Temperature (Kelvin) - * Pres = Pressure (units of if__ variable) + * Pres = Pressure (Pascal) * * Evaluate the standard state free energies at the current temperature * and pressure. Ideal gas pressure contribution is added in here. @@ -154,7 +154,7 @@ namespace VCSnonideal { for (int iph = 0; iph < NPhase; iph++) { vcs_VolPhase* vph = VPhaseList[iph]; - vph->setState_TP(m_temperature, m_pressure); + vph->setState_TP(m_temperature, m_pressurePA); vph->sendToVCSGStar(VCS_DATA_PTR(m_SSfeSpecies)); } diff --git a/Cantera/src/equil/vcs_VolPhase.cpp b/Cantera/src/equil/vcs_VolPhase.cpp index 66463aca4..ca8ea2975 100644 --- a/Cantera/src/equil/vcs_VolPhase.cpp +++ b/Cantera/src/equil/vcs_VolPhase.cpp @@ -40,6 +40,7 @@ vcs_VolPhase::vcs_VolPhase() : ElGlobalIndex(0), NVolSpecies(0), TMolesInert(0.0), + m_molarVolInert(1000.), ActivityConvention(0), Existence(0), IndexSpecialSpecies(-1), @@ -402,7 +403,7 @@ double vcs_VolPhase::G0_calc_one(int kspec, double tkelvin) { * The results are held internally within the object. * * @param TKelvin Current temperature - * @param pres Current pressure + * @param pres Current pressure (pascal) */ void vcs_VolPhase::GStar_calc(double tkelvin, double pres) { setState_TP(tkelvin, pres); @@ -431,8 +432,8 @@ void vcs_VolPhase::GStar_calc(double tkelvin, double pres) { * The kth species standard state G is returned * * @param kspec Species number (within the phase) - * @param TKelvin Current temperature - * @param pres Current pressure + * @param TKelvin Current temperature (kelvin) + * @param pres Current pressure (pascal) * * @return Gstar[kspec] returns the gibbs free energy for the * standard state of the kth species. @@ -694,7 +695,9 @@ void vcs_VolPhase::setState_TP(double temp, double pres) * The results are held internally within the object. * * @param TKelvin Current temperature - * @param pres Current pressure + * @param pres Current pressure (pascal) + * + * Calculations are in m**3/kmol */ void vcs_VolPhase::VolStar_calc(double tkelvin, double pres) { setState_TP(tkelvin, pres); @@ -722,7 +725,7 @@ void vcs_VolPhase::VolStar_calc(double tkelvin, double pres) { * * @param kspec Species number (within the phase) * @param TKelvin Current temperature - * @param pres Current pressure + * @param pres Current pressure (pascal) * * @return molar volume of the kspec species's standard * state @@ -747,8 +750,7 @@ double vcs_VolPhase::VolPM_calc() const { kglob = IndSpecies[k]; vcs_SpeciesProperties *sProp = ListSpeciesPtr[k]; VCS_SPECIES_THERMO *sTherm = sProp->SpeciesThermo; - StarMolarVol[k] = - (sTherm->VolStar_calc(kglob, Temp, Pres)); + StarMolarVol[k] = (sTherm->VolStar_calc(kglob, Temp, Pres)); } for (k = 0; k < NVolSpecies; k++) { PartialMolarVol[k] = StarMolarVol[k]; @@ -760,6 +762,16 @@ double vcs_VolPhase::VolPM_calc() const { Vol += PartialMolarVol[k] * Xmol[k]; } Vol *= TMoles; + + if (TMolesInert > 0.0) { + if (GasPhase) { + double volI = TMolesInert * 8314.47215 * Temp / Pres; + Vol += volI; + } else { + printf("unknown situation\n"); + exit(-1); + } + } } m_UpToDate_VolPM = true; return Vol; diff --git a/Cantera/src/equil/vcs_VolPhase.h b/Cantera/src/equil/vcs_VolPhase.h index 19285ebff..e6f26ebbe 100644 --- a/Cantera/src/equil/vcs_VolPhase.h +++ b/Cantera/src/equil/vcs_VolPhase.h @@ -215,6 +215,12 @@ public: //! Total moles of inert in the phase double TMolesInert; + //! molar volume of the inert species if present + /*! + * units m**3 / kg + */ + double m_molarVolInert; + //! Convention for the activity formulation /*! * 0 = molar based activities (default) @@ -285,10 +291,10 @@ public: * * Chem_Pot Pres vol moles * ---------------------------------------------------------------------- - * -1 VCS_UNITS_KCALMOL = kcal/mol atm cm**3 gmol - * 0 VCS_UNITS_UNITLESS = MU / RT -> no units atm cm**3 gmol - * 1 VCS_UNITS_KJMOL = kJ / mol atm cm**3 gmol - * 2 VCS_UNITS_KELVIN = KELVIN -> MU / R atm cm**3 gmol + * -1 VCS_UNITS_KCALMOL = kcal/gmol Pa m**3 kmol + * 0 VCS_UNITS_UNITLESS = MU / RT -> no units Pa m**3 kmol + * 1 VCS_UNITS_KJMOL = kJ / gmol Pa m**3 kmol + * 2 VCS_UNITS_KELVIN = KELVIN -> MU / R Pa m**3 kmol * 3 VCS_UNITS_MKS = Joules / Kmol (Cantera) Pa m**3 kmol * ---------------------------------------------------------------------- * @@ -358,15 +364,13 @@ public: //! Vector of the Star molar Volumes of the species. /*! - * units depends on VCS_UnitsFormat variable - * Cantera -> m3 / kmol + * units m3 / kmol */ mutable std::vector StarMolarVol; //! Vector of the Partial molar Volumes of the species. /*! - * units depends on VCS_UnitsFormat variable - * Cantera -> m3 / kmol + * units m3 / kmol */ mutable std::vector PartialMolarVol; @@ -591,6 +595,8 @@ public: * * @param TKelvin Current temperature * @param pres Current pressure + * + * Units are in m**3/kmol */ void VolStar_calc(double TKelvin, double pres); @@ -605,7 +611,7 @@ public: * @param pres Current pressure * * @return molar volume of the kspec species's standard - * state + * state (m**3/kmol) */ double VolStar_calc_one(int kglob, double TKelvin, double pres); diff --git a/Cantera/src/equil/vcs_defs.h b/Cantera/src/equil/vcs_defs.h index 1a00282f6..42046a326 100644 --- a/Cantera/src/equil/vcs_defs.h +++ b/Cantera/src/equil/vcs_defs.h @@ -191,10 +191,10 @@ namespace VCSnonideal { * @verbatim Chem_Pot Pres vol moles ------------------------------------------------- - VCS_UNITS_KCALMOL = kcal/mol atm cm**3 gmol - VCS_UNITS_UNITLESS = MU / RT -> no units atm cm**3 gmol - VCS_UNITS_KJMOL = kJ / mol atm cm**3 gmol - VCS_UNITS_KELVIN = KELVIN -> MU / R atm cm**3 gmol + VCS_UNITS_KCALMOL = kcal/mol Pa m**3 kmol + VCS_UNITS_UNITLESS = MU / RT -> no units Pa m**3 kmol + VCS_UNITS_KJMOL = kJ / mol Pa m**3 kmol + VCS_UNITS_KELVIN = KELVIN -> MU / R Pa m**3 kmol VCS_UNITS_MKS = Joules / Kmol (Cantera) Pa m**3 kmol Energy: @@ -205,10 +205,10 @@ namespace VCSnonideal { VCS_UNITS_MKS = J / kmol Pressure: (Pref and Pres) - VCS_UNITS_KCALMOL = atm - VCS_UNITS_UNITLESS = no units - VCS_UNITS_KJMOL = atm - VCS_UNITS_KELVIN = atm + VCS_UNITS_KCALMOL = Pa + VCS_UNITS_UNITLESS = Pa + VCS_UNITS_KJMOL = Pa + VCS_UNITS_KELVIN = Pa VCS_UNITS_MKS = Pa = kg / m s2 @endverbatim * @{ diff --git a/Cantera/src/equil/vcs_nondim.cpp b/Cantera/src/equil/vcs_nondim.cpp index fca355ce9..e8fe67ca2 100644 --- a/Cantera/src/equil/vcs_nondim.cpp +++ b/Cantera/src/equil/vcs_nondim.cpp @@ -113,11 +113,13 @@ void VCS_SOLVE::vcs_nondim_TP(void) { if (m_VCS_UnitsFormat == VCS_UNITS_MKS) { for (i = 0; i < m_numSpeciesTot; ++i) { if (SpeciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - m_molNumSpecies_old[i] *= 1.0E3; + //m_molNumSpecies_old[i] *= 1.0E3; + m_molNumSpecies_old[i] *= 1.0; } } for (i = 0; i < m_numElemConstraints; ++i) { - m_elemAbundancesGoal[i] *= 1.0E3; + //m_elemAbundancesGoal[i] *= 1.0E3; + m_elemAbundancesGoal[i] *= 1.0; } } } @@ -153,11 +155,13 @@ void VCS_SOLVE::vcs_redim_TP(void) if (m_VCS_UnitsFormat == VCS_UNITS_MKS) { for (i = 0; i < m_numSpeciesTot; ++i) { if (SpeciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - m_molNumSpecies_old[i] /= 1.0E3; + //m_molNumSpecies_old[i] /= 1.0E3; + m_molNumSpecies_old[i] /= 1.0; } } for (i = 0; i < m_numElemConstraints; ++i) { - m_elemAbundancesGoal[i] /= 1.0E3; + //m_elemAbundancesGoal[i] /= 1.0E3; + m_elemAbundancesGoal[i] /= 1.0; } } diff --git a/Cantera/src/equil/vcs_prep.cpp b/Cantera/src/equil/vcs_prep.cpp index 5f6c28cec..f2b4a3fd0 100644 --- a/Cantera/src/equil/vcs_prep.cpp +++ b/Cantera/src/equil/vcs_prep.cpp @@ -191,8 +191,8 @@ int VCS_SOLVE::vcs_prep_oneTime(int printLvl) } if (fabs(sum) < 1.0E-6) { modifiedSoln = true; - if (m_pressure <= 0.0) pres = 1.0; - else pres = m_pressure; + if (m_pressurePA <= 0.0) pres = 1.01325E5; + else pres = m_pressurePA; retn = vcs_evalSS_TP(0, 0, m_temperature, pres); for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) { if (SpeciesUnknownType[kspec] == VCS_SPECIES_TYPE_MOLNUM) { diff --git a/Cantera/src/equil/vcs_prob.cpp b/Cantera/src/equil/vcs_prob.cpp index a62cd4fd8..b1892580a 100644 --- a/Cantera/src/equil/vcs_prob.cpp +++ b/Cantera/src/equil/vcs_prob.cpp @@ -42,7 +42,7 @@ VCS_PROB::VCS_PROB(int nsp, int nel, int nph) : NPhase(nph), NPHASE0(0), T(298.15), - Pres(1.0), + PresPA(1.0), Vol(0.0), m_VCS_UnitsFormat(VCS_UNITS_UNITLESS), /* Set the units for the chemical potential data to be @@ -243,10 +243,8 @@ void VCS_PROB::prob_report(int print_lvl) { if (prob_type == 0) { plogf("\tSolve a constant T, P problem:\n"); plogf("\t\tT = %g K\n", T); - double pres_atm = Pres; - if (m_VCS_UnitsFormat == VCS_UNITS_MKS) { - pres_atm = Pres / 1.0133E5; - } + double pres_atm = PresPA / 1.01325E5; + plogf("\t\tPres = %g atm\n", pres_atm); } else { plogf("\tUnknown problem type\n"); @@ -279,7 +277,7 @@ void VCS_PROB::prob_report(int print_lvl) { plogf(" Information about phases\n"); plogf(" PhaseName PhaseNum SingSpec GasPhase " " EqnState NumSpec"); - plogf(" TMolesInert Tmoles\n"); + plogf(" TMolesInert TKmoles\n"); for (iphase = 0; iphase < NPhase; iphase++) { Vphase = VPhaseList[iphase]; @@ -293,10 +291,11 @@ void VCS_PROB::prob_report(int print_lvl) { } plogf("\nElemental Abundances: "); - plogf(" Target_gmol ElemType ElActive\n"); + plogf(" Target_kmol ElemType ElActive\n"); double fac = 1.0; if (m_VCS_UnitsFormat == VCS_UNITS_MKS) { - fac = 1.0E3; + //fac = 1.0E3; + fac = 1.0; } for (i = 0; i < ne; ++i) { print_space(26); plogf("%-2.2s", ElName[i].c_str()); @@ -322,7 +321,7 @@ void VCS_PROB::prob_report(int print_lvl) { for (iphase = 0; iphase < NPhase; iphase++) { Vphase = VPhaseList[iphase]; Vphase->G0_calc(T); - Vphase->GStar_calc(T, Pres); + Vphase->GStar_calc(T, PresPA); for (int kindex = 0; kindex < Vphase->NVolSpecies; kindex++) { int kglob = Vphase->IndSpecies[kindex]; plogf("%16s ", SpName[kglob].c_str()); @@ -475,7 +474,6 @@ void VCS_PROB::reportCSV(const std::string &reportFile) { exit(-1); } double Temp = T; - double pres = Pres; std::vector volPM(nspecies, 0.0); std::vector activity(nspecies, 0.0);; @@ -508,7 +506,7 @@ void VCS_PROB::reportCSV(const std::string &reportFile) { 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,"Pressure = %11.5g Pascal\n", PresPA); fprintf(FP,"Total Volume = %11.5g m**3\n", vol); fprintf(FP,"Number Basis optimizations = %d\n", m_NumBasisOptimizations); fprintf(FP,"Number VCS iterations = %d\n", m_Iterations); diff --git a/Cantera/src/equil/vcs_prob.h b/Cantera/src/equil/vcs_prob.h index e85924e51..0b6e09312 100644 --- a/Cantera/src/equil/vcs_prob.h +++ b/Cantera/src/equil/vcs_prob.h @@ -136,8 +136,9 @@ namespace VCSnonideal { //! Pressure /*! * units given by m_VCS_UnitsFormat + * -> are now PA */ - double Pres; + double PresPA; //! Volume of the entire system /*! diff --git a/Cantera/src/equil/vcs_report.cpp b/Cantera/src/equil/vcs_report.cpp index 2341b20a2..88e9bb651 100644 --- a/Cantera/src/equil/vcs_report.cpp +++ b/Cantera/src/equil/vcs_report.cpp @@ -105,23 +105,19 @@ int VCS_SOLVE::vcs_report(int iconv) * Calculate some quantities that may need updating */ vcs_tmoles(); - Vol = vcs_VolTotal(m_temperature, m_pressure, + Vol = vcs_VolTotal(m_temperature, m_pressurePA, VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(VolPM)); plogf("\t\tTemperature = %15.2g Kelvin\n", m_temperature); - std::string punits = "Atm"; - if (m_VCS_UnitsFormat == 3) { - punits = "Pa "; - } - plogf("\t\tPressure = %15.5g %3s\n", m_pressure, punits.c_str()); - plogf("\t\tVolume = %15.5g cm**3\n", Vol); + plogf("\t\tPressure = %15.5g Pa \n", m_pressurePA); + plogf("\t\tVolume = %15.5g m**3\n", Vol); /* * -------- TABLE OF SPECIES IN DECREASING MOLE NUMBERS -------------- */ plogf("\n\n"); print_line("-", 80); - plogf(" Species Equilibrium moles "); + plogf(" Species Equilibrium kmoles "); plogf("Mole Fraction ChemPot/RT SpecUnkType\n"); print_line("-", 80); for (i = 0; i < m_numComponents; ++i) { @@ -138,7 +134,7 @@ int VCS_SOLVE::vcs_report(int iconv) if (SpeciesUnknownType[l] == VCS_SPECIES_TYPE_MOLNUM) { plogf("%14.7E %14.7E %12.4E", m_molNumSpecies_old[l], m_molNumSpecies_new[l], m_feSpecies_curr[l]); - plogf(" MolNum "); + plogf(" KMolNum "); } else if (SpeciesUnknownType[l] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { plogf(" NA %14.7E %12.4E", 1.0, m_feSpecies_curr[l]); plogf(" Voltage = %14.7E", m_molNumSpecies_old[l]); @@ -162,13 +158,13 @@ int VCS_SOLVE::vcs_report(int iconv) } } if (m_numSpeciesRdc != nspecies) { - plogf("\n SPECIES WITH LESS THAN 1.0E-32 MOLES:\n\n"); + plogf("\n SPECIES WITH LESS THAN 1.0E-32 KMOLES:\n\n"); for (kspec = m_numSpeciesRdc; kspec < nspecies; ++kspec) { plogf(" %-12.12s", SpName[kspec].c_str()); plogf(" %14.7E %14.7E %12.4E", m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec], m_deltaGRxn_new[kspec]); if (SpeciesUnknownType[i] == VCS_SPECIES_TYPE_MOLNUM) { - plogf(" Mol_Num"); + plogf(" KMol_Num"); } else if (SpeciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { plogf(" Voltage"); } else { @@ -237,7 +233,7 @@ int VCS_SOLVE::vcs_report(int iconv) plogf(" %10.10s", (ElName[j]).c_str()); } plogf(" | |\n"); - plogf(" PhaseName | MolTarget |"); + plogf(" PhaseName |KMolTarget |"); for (j = 0; j < m_numElemConstraints; j++) { plogf(" %10.3g", m_elemAbundancesGoal[j]); } @@ -291,7 +287,7 @@ int VCS_SOLVE::vcs_report(int iconv) if (inertYes) plogf("\t\t(Inert species have standard free energy of zero)\n"); - plogf("\nElemental Abundances: "); + plogf("\nElemental Abundances (kmol): "); plogf(" Actual Target Type ElActive\n"); for (i = 0; i < m_numElemConstraints; ++i) { print_space(26); plogf("%-2.2s", (ElName[i]).c_str()); @@ -310,7 +306,7 @@ int VCS_SOLVE::vcs_report(int iconv) plogf("\t\t(RT = %g ", rt); vcs_printChemPotUnits(m_VCS_UnitsFormat); plogf(")\n"); - plogf(" Name TMoles StandStateChemPot " + plogf(" Name TKMoles StandStateChemPot " " ln(AC) ln(X_i) | F z_i phi | ChemPot | (-lnMnaught)\n"); print_line("-", 115); for (i = 0; i < nspecies; ++i) { diff --git a/Cantera/src/equil/vcs_solve.cpp b/Cantera/src/equil/vcs_solve.cpp index 0168c5f4c..09eee6c83 100644 --- a/Cantera/src/equil/vcs_solve.cpp +++ b/Cantera/src/equil/vcs_solve.cpp @@ -44,7 +44,7 @@ namespace VCSnonideal { iest(0), TMoles(0.0), m_temperature(0.0), - m_pressure(0.0), + m_pressurePA(0.0), tolmaj(0.0), tolmin(0.0), tolmaj2(0.0), @@ -389,7 +389,7 @@ namespace VCSnonideal { * calculate the residual and Jacobian) */ - iconv = vcs_TP(ipr, ip1, maxit, vprob->T, vprob->Pres); + iconv = vcs_TP(ipr, ip1, maxit, vprob->T, vprob->PresPA); /* @@ -582,8 +582,8 @@ namespace VCSnonideal { */ if (pub->T > 0.0) m_temperature = pub->T; else m_temperature = 293.15; - if (pub->Pres > 0.0) m_pressure = pub->Pres; - else m_pressure = 1.0; + if (pub->PresPA > 0.0) m_pressurePA = pub->PresPA; + else m_pressurePA = Cantera::OneAtm; /* * TPhInertMoles[] -> must be copied over here */ @@ -773,7 +773,7 @@ namespace VCSnonideal { bool status_change = false; m_temperature = pub->T; - m_pressure = pub->Pres; + m_pressurePA = pub->PresPA; m_VCS_UnitsFormat = pub->m_VCS_UnitsFormat; iest = pub->iest; @@ -902,7 +902,7 @@ namespace VCSnonideal { int k1 = 0; vcs_tmoles(); - Vol = vcs_VolTotal(m_temperature, m_pressure, + Vol = vcs_VolTotal(m_temperature, m_pressurePA, VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(VolPM)); for (i = 0; i < m_numSpeciesTot; ++i) { @@ -930,7 +930,7 @@ namespace VCSnonideal { } pub->T = m_temperature; - pub->Pres = m_pressure; + pub->PresPA = m_pressurePA; pub->Vol = Vol; int kT = 0; for (int iph = 0; iph < pub->NPhase; iph++) { diff --git a/Cantera/src/equil/vcs_solve.h b/Cantera/src/equil/vcs_solve.h index d27595476..878c977ab 100644 --- a/Cantera/src/equil/vcs_solve.h +++ b/Cantera/src/equil/vcs_solve.h @@ -621,7 +621,7 @@ public: * Note other constraint conditions are added to this vector. * This is input from the input file and * is considered a constant from thereon. - * units = gmoles + * units = kmoles */ std::vector m_elemAbundances; @@ -632,7 +632,7 @@ public: * Note other constraint conditions are added to this vector. * This is input from the input file and * is considered a constant from thereon. - * units = gmoles + * units = kmoles */ std::vector m_elemAbundancesGoal; @@ -641,7 +641,7 @@ public: * -> Don't use this except for scaling * purposes only */ - //! total gmols of species in each phase + //! total kmols of species in each phase /*! * This contains the total number of moles of species in each phase * @@ -649,7 +649,7 @@ public: */ std::vector m_tPhaseMoles_old; - //! total gmols of species in each phase in the tentative soln vector + //! total kmols of species in each phase in the tentative soln vector /*! * This contains the total number of moles of species in each phase * in the tentative solution vector @@ -681,12 +681,13 @@ public: * 1: atm * 2: atm * 3: Pa + * Units being changed to Pa */ - double m_pressure; + double m_pressurePA; //! Total kmoles of inert to add to each phase /*! - * TPhInertMoles[iph] = Total gmoles of inert to add to each phase + * TPhInertMoles[iph] = Total kmoles of inert to add to each phase * length = number of phases */ std::vector TPhInertMoles; @@ -864,7 +865,7 @@ public: //! Molecular weight of each species /*! - * units = gm/gmol + * units = kg/kmol * length = number of species */ std::vector WtSpecies; @@ -891,7 +892,7 @@ public: */ int UseActCoeffJac; - double Vol; /* Vol = Volume (cm^3) */ + double Vol; /* Vol = Volume (m^3) */ //! Partialm molar volumes of the species /*! @@ -928,10 +929,10 @@ public: * 2: Kelvin * 3: J / kmol * and pressure data: - * -1: atm - * 0: atm - * 1: atm - * 2: atm + * -1: Pa + * 0: Pa + * 1: Pa + * 2: pa * 3: Pa */ int m_VCS_UnitsFormat; diff --git a/Cantera/src/equil/vcs_solve_TP.cpp b/Cantera/src/equil/vcs_solve_TP.cpp index 4332f5dc9..caa209a76 100644 --- a/Cantera/src/equil/vcs_solve_TP.cpp +++ b/Cantera/src/equil/vcs_solve_TP.cpp @@ -237,11 +237,11 @@ namespace VCSnonideal { plogf("%10d PHASE2 SPECIES%8d SINGLE SPECIES PHASES\n\n", numSpecliquid, m_numSpeciesTot - (VPhaseList[0])->NVolSpecies - numSpecliquid); - string punits = "atm"; - if (m_VCS_UnitsFormat == 3) { - punits = "Pa "; - } - plogf(" PRESSURE%22.8g %3s\n", m_pressure, punits.c_str()); + //string punits = "atm"; + //if (m_VCS_UnitsFormat == 3) { + // punits = "Pa "; + //} + plogf(" PRESSURE%22.8g %3s\n", m_pressurePA, "Pa "); plogf(" TEMPERATURE%19.3f K\n", m_temperature); Vphase = VPhaseList[0]; if (Vphase->NVolSpecies > 0) { @@ -532,7 +532,7 @@ namespace VCSnonideal { if (iti == 0) plogf("- Full Calculation:\n"); else plogf("- Major Components Calculation:\n"); plogf(" --- Species IC "); - plogf(" Moles Tent_Moles Rxn_Adj | Comment \n"); + plogf(" KMoles Tent_KMoles Rxn_Adj | Comment \n"); } #endif @@ -785,7 +785,7 @@ namespace VCSnonideal { */ if (m_molNumSpecies_new[kspec] <= 0.0) { #ifdef DEBUG_MODE - sprintf(ANOTE, "initial nonpos moles= %11.3E", + sprintf(ANOTE, "initial nonpos kmoles= %11.3E", m_molNumSpecies_new[kspec]); #endif /* ************************************************* */ @@ -801,7 +801,7 @@ namespace VCSnonideal { if (! (SSPhase[kspec])) { /* * Section for multispecies phases: - * - Cut reaction adjustment for positive moles of + * - Cut reaction adjustment for positive kmoles of * major species in multispecies phases. * Decrease its concentration by a factor of 10. */ @@ -811,7 +811,7 @@ namespace VCSnonideal { /* * Change major to minor if the current species * has a mole number that is less than 1/100 of the - * total moles in the problem. + * total kmoles in the problem. * However, it also has to be a small species within its * own phase as well. * we can't call vcs_species_type() because the phase moles @@ -954,7 +954,7 @@ namespace VCSnonideal { } /* End of Loop on ic[irxn] -> the type of species */ /***********************************************************************/ - /****** CALCULATE MOLE NUMBER CHANGE FOR THE COMPONENT BASIS ***********/ + /****** CALCULATE KMOLE NUMBER CHANGE FOR THE COMPONENT BASIS **********/ /***********************************************************************/ if (dx != 0.0 && (SpeciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE)) { /* @@ -1143,7 +1143,7 @@ namespace VCSnonideal { plogf(" -----------------------------------------------------\n"); plogf(" --- FORCER SUBROUTINE changed the solution:\n"); plogf(" --- SPECIES Status INIT MOLES TENT_MOLES"); - plogf(" FINAL MOLES INIT_DEL_G/RT TENT_DEL_G/RT FINAL_DELTA_G/RT\n"); + plogf(" FINAL KMOLES INIT_DEL_G/RT TENT_DEL_G/RT FINAL_DELTA_G/RT\n"); for (i = 0; i < m_numComponents; ++i) { plogf(" --- %-12.12s", SpName[i].c_str()); plogf(" %14.6E %14.6E %14.6E\n", m_molNumSpecies_old[i], m_molNumSpecies_old[i] + m_deltaMolNumSpecies[i], m_molNumSpecies_new[i]); @@ -1160,11 +1160,11 @@ namespace VCSnonideal { plogf("Norms of Delta G():%14.6E%14.6E\n", l2normdg(VCS_DATA_PTR(m_deltaGRxn_old)), l2normdg(VCS_DATA_PTR(m_deltaGRxn_new))); - plogf(" Total moles of gas = %15.7E\n", m_tPhaseMoles_old[0]); + plogf(" Total kmoles of gas = %15.7E\n", m_tPhaseMoles_old[0]); if ((NPhase > 1) && (! (VPhaseList[1])->SingleSpecies)) { - plogf(" Total moles of liquid = %15.7E\n", m_tPhaseMoles_old[1]); + plogf(" Total kmoles of liquid = %15.7E\n", m_tPhaseMoles_old[1]); } else { - plogf(" Total moles of liquid = %15.7E\n", 0.0); + plogf(" Total kmoles of liquid = %15.7E\n", 0.0); } plogf(" Total New Dimensionless Gibbs Free Energy = %20.13E\n", vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_new), VCS_DATA_PTR(m_feSpecies_curr), @@ -1185,7 +1185,7 @@ namespace VCSnonideal { } else { plogf(" (only major species):\n"); } - plogf(" --- Species Status Initial_Moles Final_Moles Initial_Mu/RT"); + plogf(" --- Species Status Initial_KMoles Final_KMoles Initial_Mu/RT"); plogf(" Mu/RT Init_Del_G/RT Delta_G/RT\n"); for (i = 0; i < m_numComponents; ++i) { plogf(" --- %-12.12s", SpName[i].c_str()); plogf(" "); @@ -1214,7 +1214,7 @@ namespace VCSnonideal { l2normdg(VCS_DATA_PTR(m_deltaGRxn_new))); plogendl(); - plogf(" --- Phase_Name Moles(after update)\n"); + plogf(" --- Phase_Name KMoles(after update)\n"); plogf(" --- "); vcs_print_line("-", 50); for (iph = 0; iph < NPhase; iph++) { Vphase = VPhaseList[iph]; @@ -2129,7 +2129,7 @@ namespace VCSnonideal { * * Change the concentration of a species by delta moles. * Make sure to conserve - * elements and keep track of the total moles in all phases. + * elements and keep track of the total kmoles in all phases. * * return: * 1: succeeded @@ -2831,7 +2831,7 @@ namespace VCSnonideal { plogf(" "); for (j = 0; j < 82; j++) plogf("-"); plogf("\n"); plogf(" --- Subroutine vcs_RxnStepSizes called - Details:\n"); plogf(" "); for (j = 0; j < 82; j++) plogf("-"); plogf("\n"); - plogf(" --- Species Moles Rxn_Adjustment DeltaG" + plogf(" --- Species KMoles Rxn_Adjustment DeltaG" " | Comment\n"); } #endif diff --git a/Cantera/src/equil/vcs_species_thermo.cpp b/Cantera/src/equil/vcs_species_thermo.cpp index e1e463bc1..c3baace58 100644 --- a/Cantera/src/equil/vcs_species_thermo.cpp +++ b/Cantera/src/equil/vcs_species_thermo.cpp @@ -43,7 +43,7 @@ VCS_SPECIES_THERMO::VCS_SPECIES_THERMO(int indexPhase, SS0_H0(0.0), SS0_S0(0.0), SS0_Cp0(0.0), - SS0_Pref(1.0), + SS0_Pref(1.01325E5), SS0_Params(0), SSStar_Model(VCS_SSSTAR_CONSTANT), SSStar_Params(0), @@ -55,15 +55,7 @@ VCS_SPECIES_THERMO::VCS_SPECIES_THERMO(int indexPhase, UseCanteraCalls(false), m_VCS_UnitsFormat(VCS_UNITS_UNITLESS) { - /* - * Set up the numerical value for P_reference, based on the current - * global units choice. - */ - if (m_VCS_UnitsFormat == VCS_UNITS_MKS) { - SS0_Pref = 1.0133E5; - } else { - SS0_Pref = 1.0; - } + SS0_Pref = 1.01325E5; } @@ -224,11 +216,11 @@ double VCS_SPECIES_THERMO::GStar_R_calc(int kglob, double TKelvin, * Input * * Output - * return value = standard state volume in cm**3 per mol. - * (VCS_UNITS_MKS) m**3 / kmol + * return value = standard state volume in m**3 per kmol. + * (VCS_UNITS_MKS) */ double VCS_SPECIES_THERMO:: -VolStar_calc(int kglob, double TKelvin, double pres) +VolStar_calc(int kglob, double TKelvin, double presPA) { char yo[] = "VCS_SPECIES_THERMO::VStar_calc "; double vol, T; @@ -237,18 +229,15 @@ VolStar_calc(int kglob, double TKelvin, double pres) if (UseCanteraCalls) { AssertThrowVCS(m_VCS_UnitsFormat == VCS_UNITS_MKS, "Possible inconsistency"); int kspec = IndexSpeciesPhase; - vol = OwningPhase->VolStar_calc_one(kspec, TKelvin, pres); + vol = OwningPhase->VolStar_calc_one(kspec, TKelvin, presPA); } else { switch(SSStar_Vol_Model) { case VCS_SSVOL_CONSTANT: vol = SSStar_Vol0; break; case VCS_SSVOL_IDEALGAS: - if (m_VCS_UnitsFormat == VCS_UNITS_MKS) { - vol = 8.31451E3 * T / pres; - } else { - vol= 83.14510 / 1.01325 * T / pres; - } + // R J/kmol/K (2006 CODATA value) + vol= 8314.47215 * T / presPA; break; default: plogf("%sERROR: unknown SSVol model\n", yo); diff --git a/Cantera/src/equil/vcs_species_thermo.h b/Cantera/src/equil/vcs_species_thermo.h index e7e68e9ca..268b9731f 100644 --- a/Cantera/src/equil/vcs_species_thermo.h +++ b/Cantera/src/equil/vcs_species_thermo.h @@ -112,6 +112,7 @@ public: /** * Value of the pressure for the reference state. + * defaults to 1.01325E5 = 1 atm */ double SS0_Pref; /**