diff --git a/test_problems/SConscript b/test_problems/SConscript index d4f926301..0b6dde927 100644 --- a/test_problems/SConscript +++ b/test_problems/SConscript @@ -267,6 +267,8 @@ CompileAndTest('rankine_democxx', 'rankine_democxx', 'rankine', 'output_blessed. artifacts=['liquidvapor.xml']) CompileAndTest('silane_equil', 'silane_equil', 'silane_equi', 'output_blessed.txt') # spectroscopy is incomplete +CompileAndTest('stoichSolidKinetics', 'stoichSolidKinetics', + 'stoichSolidKinetics', 'output_blessed.txt') CompileAndTest('surfkin', 'surfkin', 'surfdemo', 'output_blessed.txt') CompileAndTest('surfSolver', 'surfSolverTest', 'surfaceSolver', None, arguments='haca2.xml', diff --git a/test_problems/stoichSolidKinetics/ReactionSurf.xml b/test_problems/stoichSolidKinetics/ReactionSurf.xml new file mode 100644 index 000000000..6458c8f73 --- /dev/null +++ b/test_problems/stoichSolidKinetics/ReactionSurf.xml @@ -0,0 +1,71 @@ + + + + + + + + Ca C O Fe N + empty_site + + + 1500.0 + empty_site:1.0 + + + 3e-09 + + + + + + CaCO3(S) air CaO(S) Fe3O4(S) FeO(S) Fe(S) + + + + + + + + + + + + 298.15 + 0.0 + 0.0 + 0.0 + + + + + + + + + + + CaCO3(S) [=] CaO(S) + CO2 + + + 1.0E-5 + 0.0 + 20.000000 + + + CaCO3(S):1.0 + CaO(S):1.0 CO2:1.0 + + + + + + diff --git a/test_problems/stoichSolidKinetics/gas.xml b/test_problems/stoichSolidKinetics/gas.xml new file mode 100644 index 000000000..dc2e7f15d --- /dev/null +++ b/test_problems/stoichSolidKinetics/gas.xml @@ -0,0 +1,252 @@ + + + + + + + + O H C Fe Ca N + + + O2 H2 H CH2O CHO OH CO2 H2O CH4 CO N2 + + + + 1500.0 + 101325.0 + + + + + + + + + + + + + + O:2 + + + 1500.0 + 0.0 + 0.0 + 0.0 + + + + linear + 107.400 + 3.460 + 0.000 + 1.600 + 3.800 + + + + + H:2 + + + 1500.0 + 0.0 + 0.0 + 0.0 + + + + linear + 107.400 + 3.460 + 0.000 + 1.600 + 3.800 + + + + + H:1 + + + 1500.0 + 1.0E6 + 0.0 + 0.0 + + + + linear + 107.400 + 3.460 + 0.000 + 1.600 + 3.800 + + + + + C:1 H:2 O:1 + + + 1500.0 + -8.611E4 + 0.0 + 0.0 + + + + linear + 107.400 + 3.460 + 0.000 + 1.600 + 3.800 + + + + + C:1 H:1 O:1 + + + 1500.0 + -6.255E4 + 0.0 + 0.0 + + + + linear + 107.400 + 3.460 + 0.000 + 1.600 + 3.800 + + + + + H:1 O:1 + + + 1500.0 + 2.259E4 + 0.0 + 0.0 + + + + linear + 107.400 + 3.460 + 0.000 + 1.600 + 3.800 + + + + + C:1 O:2 + + + 1500.0 + -3.9597E5 + 0.0 + 0.0 + + + + linear + 107.400 + 3.460 + 0.000 + 1.600 + 3.800 + + + + + H:2 O:1 + + + 1500.0 + -1.8987E5 + 0.0 + 0.0 + + + + linear + 107.400 + 3.460 + 0.000 + 1.600 + 3.800 + + + + + C:1 H:4 + + + 1500.0 + 2.485E4 + 0.0 + 0.0 + + + + linear + 107.400 + 3.460 + 0.000 + 1.600 + 3.800 + + + + + C:1 O:1 + + + 1500.0 + -2.0464E5 + 0.0 + 0.0 + + + + linear + 107.400 + 3.460 + 0.000 + 1.600 + 3.800 + + + + + + + N:2 + + + 1500.0 + 0.0 + 0.0 + 0.0 + + + + linear + 97.530 + 3.620 + 0.000 + 1.760 + 4.000 + + + + + diff --git a/test_problems/stoichSolidKinetics/output_blessed.txt b/test_problems/stoichSolidKinetics/output_blessed.txt new file mode 100644 index 000000000..b07993814 --- /dev/null +++ b/test_problems/stoichSolidKinetics/output_blessed.txt @@ -0,0 +1,96 @@ +*** StoichSolidKinetics Test *** +Tests for the proper behavior of heterogeneous reactions +when phases may or may not exist: + CaCO3(s) = CO2(g) + CaO(s) + + deltaGSS = 1.729E+07 + deltaG = 3.908E+06 + mu_CO2(g) = -4.094E+08 + mu_cao(s) = -5.292E+08 + mu_caco3 = -9.424E+08 + deltaG_calc = 3.908E+06 + act_CO2 = 2.000E-01 + act_cao(s) = 1.000E+00 + act_caco3(s) = 1.000E+00 +*** Base problem assuming that all phases exist: +ROP_net = -5.414E-07 +ROP_forward = 9.023E-07 +ROP_reverse = 1.444E-06 + kfwd = 9.023E-07 + krev = 5.923E-04 +*** Setting CaO(S) phase to nonexistent: +ROP_net = 0.000E+00 +ROP_forward = 9.023E-07 +ROP_reverse = 9.023E-07 + kfwd = 9.023E-07 + krev = 5.923E-04 +*** Setting CaCO3(S) phase to nonexistent: +ROP_net = 0.000E+00 +ROP_forward = 0.000E+00 +ROP_reverse = 0.000E+00 + kfwd = 9.023E-07 + krev = 5.923E-04 +*** OK Setting CaO(S) phase to existent, CaCO3 nonexistent: +ROP_net = -5.414E-07 +ROP_forward = 9.023E-07 +ROP_reverse = 1.444E-06 + kfwd = 9.023E-07 + krev = 5.923E-04 +*** Setting Gas phase to nonexistent, CaCO3 nonexistent: +ROP_net = 0.000E+00 +ROP_forward = 0.000E+00 +ROP_reverse = 0.000E+00 + kfwd = 9.023E-07 + krev = 5.923E-04 +*** Setting to all phases existing: +ROP_net = -5.414E-07 +ROP_forward = 9.023E-07 +ROP_reverse = 1.444E-06 + kfwd = 9.023E-07 + krev = 5.923E-04 +*** Setting so that forward rate if faster: + deltaGSS = 1.729E+07 + deltaG = -3.438E+07 + mu_CO2(g) = -4.476E+08 + mu_cao(s) = -5.292E+08 + mu_caco3 = -9.424E+08 + deltaG_calc = -3.438E+07 + act_CO2 = 2.000E-03 + act_cao(s) = 1.000E+00 + act_caco3(s) = 1.000E+00 +*** Base problem assuming that all phases exist: +ROP_net = 8.878E-07 +ROP_forward = 9.023E-07 +ROP_reverse = 1.444E-08 + kfwd = 9.023E-07 + krev = 5.923E-04 +*** Setting CaO(S) phase to nonexistent: +ROP_net = 8.878E-07 +ROP_forward = 9.023E-07 +ROP_reverse = 1.444E-08 + kfwd = 9.023E-07 + krev = 5.923E-04 +*** Setting CaCO3(S) phase to nonexistent: +ROP_net = 0.000E+00 +ROP_forward = 0.000E+00 +ROP_reverse = 0.000E+00 + kfwd = 9.023E-07 + krev = 5.923E-04 +*** OK Setting CaO(S) phase to existent, CaCO3 nonexistent: +ROP_net = 0.000E+00 +ROP_forward = 1.444E-08 +ROP_reverse = 1.444E-08 + kfwd = 9.023E-07 + krev = 5.923E-04 +*** Setting Gas phase to nonexistent, CaCO3 nonexistent: +ROP_net = 0.000E+00 +ROP_forward = 0.000E+00 +ROP_reverse = 0.000E+00 + kfwd = 9.023E-07 + krev = 5.923E-04 +*** Setting to all phases existing: +ROP_net = 8.878E-07 +ROP_forward = 9.023E-07 +ROP_reverse = 1.444E-08 + kfwd = 9.023E-07 + krev = 5.923E-04 diff --git a/test_problems/stoichSolidKinetics/solidPhases.xml b/test_problems/stoichSolidKinetics/solidPhases.xml new file mode 100644 index 000000000..3adf45aba --- /dev/null +++ b/test_problems/stoichSolidKinetics/solidPhases.xml @@ -0,0 +1,155 @@ + + + + + + + + O H C Fe Ca N + + CaO(S) + + 3.52 + + + + + + + + O H C Fe Ca N + + CaCO3(S) + + 3.52 + + + + + + + + O H C Fe Ca N + + C(d) + + 3.52 + + + + + + + + O H C Fe Ca N + + Fe3O4(S) + + 3.52 + + + + + + + + O H C Fe Ca N + + FeO_Solid + + 3.52 + + + + + + + + O H C Fe Ca N + + Fe_Solid + + 3.52 + + + + + + + + + + O:1 Ca:1 + + + 298.14999999999998 + -5.2919E5 + 0.0 + 0.0 + + + + + + Ca:1 C:1 O:3 + + + 298.14999999999998 + -9.4245E5 + 0.0 + 0.0 + + + + + + C:1 + + + 298.14999999999998 + 0.0 + 0.0 + 0.0 + + + + + + O:4 Fe:3 + + + 298.14999999999998 + -7.6266E5 + 0.0 + 0.0 + + + + + + Fe:1 O:1 + + + 298.14999999999998 + -1.9393E5 + 0.0 + 0.0 + + + + + + Fe:1 + + + 298.14999999999998 + 0.0 + 0.0 + 0.0 + + + + + + + diff --git a/test_problems/stoichSolidKinetics/stoichSolidKinetics.cpp b/test_problems/stoichSolidKinetics/stoichSolidKinetics.cpp new file mode 100644 index 000000000..11366922b --- /dev/null +++ b/test_problems/stoichSolidKinetics/stoichSolidKinetics.cpp @@ -0,0 +1,180 @@ +#include "cantera/electrolyteThermo.h" +#include "cantera/kinetics.h" +#include "cantera/thermo/MargulesVPSSTP.h" +#include "cantera/thermo/IonsFromNeutralVPSSTP.h" +#include "cantera/IdealGasMix.h" +#include "cantera/base/PrintCtrl.h" + +#include + +using namespace Cantera; +using namespace std; + +Cantera::PrintCtrl pc; + +void printValue(const std::string& label, double value) +{ + std::cout << label; + pc.pr_de(value, 4, 13); + std::cout << std::endl; +} + +void printRates(InterfaceKinetics& iKin) +{ + vector_fp work(iKin.nReactions(), 0.0); + iKin.getNetRatesOfProgress(&work[0]); + printValue("ROP_net = ", work[0]); + + iKin.getFwdRatesOfProgress(&work[0]); + printValue("ROP_forward = ", work[0]); + + iKin.getRevRatesOfProgress(&work[0]); + printValue("ROP_reverse = ", work[0]); + + iKin.getFwdRateConstants(&work[0]); + printValue(" kfwd = ", work[0]); + + iKin.getRevRateConstants(&work[0]); + printValue(" krev = ", work[0]); +} + +void testProblem() +{ + XML_Node xc; + ctml::get_CTML_Tree(&xc, "ReactionSurf.xml"); + XML_Node* xg = xc.findNameID("phase", "reaction_surface"); + if (!xg) { + throw CanteraError("couldn't find file", ""); + } + std::auto_ptr surfTP(newPhase(*xg)); + std::auto_ptr gasTP(newPhase("gas.xml", "")); + + std::auto_ptr cao_s(newPhase("solidPhases.xml", "CaO(S)")); + std::auto_ptr caco3_s(newPhase("solidPhases.xml", "CaCO3(S)")); + std::auto_ptr c_s(newPhase("solidPhases.xml", "C(S)")); + std::auto_ptr fe3o4_s(newPhase("solidPhases.xml", "Fe3O4(S)")); + std::auto_ptr feo_s(newPhase("solidPhases.xml", "FeO(S)")); + std::auto_ptr fe_s(newPhase("solidPhases.xml", "Fe(S)")); + + vector phaseList; + phaseList.push_back(gasTP.get()); + phaseList.push_back(cao_s.get()); + phaseList.push_back(caco3_s.get()); + phaseList.push_back(c_s.get()); + phaseList.push_back(fe3o4_s.get()); + phaseList.push_back(feo_s.get()); + phaseList.push_back(fe_s.get()); + phaseList.push_back(surfTP.get()); + + InterfaceKinetics iKin; + importKinetics(*xg, phaseList, &iKin); + + vector_fp mll(gasTP->nSpecies(), 0.0); + int igco2 = gasTP->speciesIndex("CO2"); + int igo2 = gasTP->speciesIndex("O2"); + int ign2 = gasTP->speciesIndex("N2"); + + mll[igco2] = 0.2; + mll[igo2] = 0.1; + mll[ign2] = 0.7; + + // Set the bath gas of 1000 K and 1 atm + double Temp = 1000.; + gasTP->setState_TPX(Temp, OneAtm, &mll[0]); + cao_s->setState_TP(Temp, OneAtm); + caco3_s->setState_TP(Temp, OneAtm); + c_s->setState_TP(Temp, OneAtm); + fe3o4_s->setState_TP(Temp, OneAtm); + surfTP->setState_TP(Temp, OneAtm); + + vector_fp work(gasTP->nSpecies(), 0.0); + + cout << "*** StoichSolidKinetics Test ***" << endl; + + cout << "Tests for the proper behavior of heterogeneous reactions\n" + << "when phases may or may not exist:\n" + << " CaCO3(s) = CO2(g) + CaO(s) \n" << endl; + + for (int ktrials = 0; ktrials < 2; ktrials++) { + iKin.getDeltaSSGibbs(&work[0]); + printValue(" deltaGSS = ", work[0]); + + iKin.getDeltaGibbs(&work[0]); + printValue(" deltaG = ", work[0]); + + gasTP->getChemPotentials(&work[0]); + double mu_CO2 = work[igco2]; + printValue(" mu_CO2(g) = ", mu_CO2); + + cao_s->getGibbs_RT(&work[0]); + double mu_cao = work[0] * GasConstant * Temp; + printValue(" mu_cao(s) = ", mu_cao); + + caco3_s->getChemPotentials(&work[0]); + double mu_caco3 = work[0]; + printValue(" mu_caco3 = ", mu_caco3); + + double deltaG_calc = mu_CO2 + mu_cao - mu_caco3; + printValue(" deltaG_calc = ", deltaG_calc); + + gasTP->getActivities(&work[0]); + double act_CO2 = work[igco2]; + printValue(" act_CO2 = ", act_CO2); + + cao_s->getActivities(&work[0]); + printValue(" act_cao(s) = ", work[0]); + + caco3_s->getActivities(&work[0]); + printValue(" act_caco3(s) = ", work[0]); + + cout << "*** Base problem assuming that all phases exist:" << endl; + printRates(iKin); + + cout << "*** Setting CaO(S) phase to nonexistent:" << endl; + int ip_cao = iKin.phaseIndex("CaO(S)"); + iKin.setPhaseExistence(ip_cao, false); + iKin.setPhaseStability(ip_cao, true); + printRates(iKin); + + cout << "*** Setting CaCO3(S) phase to nonexistent:" << endl; + int ip_caco3 = iKin.phaseIndex("CaCO3(S)"); + iKin.setPhaseExistence(ip_caco3, false); + iKin.setPhaseStability(ip_caco3, true); + printRates(iKin); + + cout << "*** OK Setting CaO(S) phase to existent, CaCO3 nonexistent:" << endl; + iKin.setPhaseExistence(ip_cao, true); + printRates(iKin); + + cout << "*** Setting Gas phase to nonexistent, CaCO3 nonexistent:" << endl; + int ip_gas = iKin.phaseIndex("air"); + iKin.setPhaseExistence(ip_gas, false); + iKin.setPhaseStability(ip_gas, true); + printRates(iKin); + + cout << "*** Setting to all phases existing:" << endl; + iKin.setPhaseExistence(ip_gas, true); + iKin.setPhaseExistence(ip_cao, true); + iKin.setPhaseExistence(ip_caco3, true); + printRates(iKin); + + if (ktrials == 0) { + cout << "*** Setting so that forward rate if faster:" << endl; + mll[igco2] = 0.002; + mll[igo2] = 0.1; + mll[ign2] = 0.898; + gasTP->setState_TPX(Temp, OneAtm, &mll[0]); + } + } +} + +int main(int argc, char** argv) +{ + try { + testProblem(); + return 0; + } catch (CanteraError& err) { + std::cout << err.what() << std::endl; + return 0; + } +}