Cantera  3.1.0a1
Loading...
Searching...
No Matches
vcs_MultiPhaseEquil.cpp
Go to the documentation of this file.
1/**
2 * @file vcs_MultiPhaseEquil.cpp
3 * Driver routine for the VCSnonideal equilibrium solver package
4 */
5
6// This file is part of Cantera. See License.txt in the top-level directory or
7// at https://cantera.org/license.txt for license and copyright information.
8
17
18#include <cstdio>
19
20namespace Cantera
21{
23 m_mix(mix),
24 m_printLvl(printLvl),
25 m_vsolve(mix, printLvl)
26{
27}
28
29int vcs_MultiPhaseEquil::equilibrate_TV(int XY, double xtarget, int estimateEquil,
30 int printLvl, double err,
31 int maxsteps, int loglevel)
32{
33 double Vtarget = m_mix->volume();
34 if ((XY != TV) && (XY != HV) && (XY != UV) && (XY != SV)) {
35 throw CanteraError("vcs_MultiPhaseEquil::equilibrate_TV",
36 "Wrong XY flag: {}", XY);
37 }
38 int maxiter = 100;
39 int iSuccess = 0;
40 if (XY == TV) {
41 m_mix->setTemperature(xtarget);
42 }
43 int strt = estimateEquil;
44 double P1 = 0.0;
45 double V1 = 0.0;
46 double V2 = 0.0;
47 double P2 = 0.0;
48 double Tlow = 0.5 * m_mix->minTemp();
49 double Thigh = 2.0 * m_mix->maxTemp();
50 int printLvlSub = std::max(0, printLvl - 1);
51 for (int n = 0; n < maxiter; n++) {
52 double Pnow = m_mix->pressure();
53
54 switch (XY) {
55 case TV:
56 iSuccess = equilibrate_TP(strt, printLvlSub, err, maxsteps, loglevel);
57 break;
58 case HV:
59 iSuccess = equilibrate_HP(xtarget, HP, Tlow, Thigh, strt,
60 printLvlSub, err, maxsteps, loglevel);
61 break;
62 case UV:
63 iSuccess = equilibrate_HP(xtarget, UP, Tlow, Thigh, strt,
64 printLvlSub, err, maxsteps, loglevel);
65 break;
66 case SV:
67 iSuccess = equilibrate_SP(xtarget, Tlow, Thigh, strt,
68 printLvlSub, err, maxsteps, loglevel);
69 break;
70 default:
71 break;
72 }
73 strt = false;
74 double Vnow = m_mix->volume();
75 if (n == 0) {
76 V2 = Vnow;
77 P2 = Pnow;
78 } else if (n == 1) {
79 V1 = Vnow;
80 P1 = Pnow;
81 } else {
82 P2 = P1;
83 V2 = V1;
84 P1 = Pnow;
85 V1 = Vnow;
86 }
87
88 double Verr = fabs((Vtarget - Vnow)/Vtarget);
89 if (Verr < err) {
90 return iSuccess;
91 }
92 double Pnew;
93 // find dV/dP
94 if (n > 1) {
95 double dVdP = (V2 - V1) / (P2 - P1);
96 if (dVdP == 0.0) {
97 throw CanteraError("vcs_MultiPhase::equilibrate_TV",
98 "dVdP == 0.0");
99 } else {
100 Pnew = Pnow + (Vtarget - Vnow) / dVdP;
101 if (Pnew < 0.2 * Pnow) {
102 Pnew = 0.2 * Pnow;
103 }
104 if (Pnew > 3.0 * Pnow) {
105 Pnew = 3.0 * Pnow;
106 }
107 }
108 } else {
109 m_mix->setPressure(Pnow*1.01);
110 double dVdP = (m_mix->volume() - Vnow)/(0.01*Pnow);
111 Pnew = Pnow + 0.5*(Vtarget - Vnow)/dVdP;
112 if (Pnew < 0.5* Pnow) {
113 Pnew = 0.5 * Pnow;
114 }
115 if (Pnew > 1.7 * Pnow) {
116 Pnew = 1.7 * Pnow;
117 }
118 }
119 m_mix->setPressure(Pnew);
120 }
121 throw CanteraError("vcs_MultiPhase::equilibrate_TV",
122 "No convergence for V");
123}
124
125int vcs_MultiPhaseEquil::equilibrate_HP(double Htarget, int XY, double Tlow,
126 double Thigh, int estimateEquil, int printLvl, double err, int maxsteps,
127 int loglevel)
128{
129 int maxiter = 100;
130 int iSuccess;
131 if (XY != HP && XY != UP) {
132 throw CanteraError("vcs_MultiPhaseEquil::equilibrate_HP",
133 "Wrong XP", XY);
134 }
135 int strt = estimateEquil;
136
137 // Lower bound on T. This will change as we progress in the calculation
138 if (Tlow <= 0.0) {
139 Tlow = 0.5 * m_mix->minTemp();
140 }
141 // Upper bound on T. This will change as we progress in the calculation
142 if (Thigh <= 0.0 || Thigh > 1.0E6) {
143 Thigh = 2.0 * m_mix->maxTemp();
144 }
145
146 double cpb = 1.0;
147 double Hlow = Undef;
148 double Hhigh = Undef;
149 double Tnow = m_mix->temperature();
150 int printLvlSub = std::max(printLvl - 1, 0);
151
152 for (int n = 0; n < maxiter; n++) {
153 // start with a loose error tolerance, but tighten it as we get
154 // close to the final temperature
155 try {
156 Tnow = m_mix->temperature();
157 iSuccess = equilibrate_TP(strt, printLvlSub, err, maxsteps, loglevel);
158 strt = 0;
159 double Hnow = (XY == UP) ? m_mix->IntEnergy() : m_mix->enthalpy();
160 double pmoles[10];
161 pmoles[0] = m_mix->phaseMoles(0);
162 double Tmoles = pmoles[0];
163 double HperMole = Hnow/Tmoles;
164 if (printLvl > 0) {
165 plogf("T = %g, Hnow = %g ,Tmoles = %g, HperMole = %g\n",
166 Tnow, Hnow, Tmoles, HperMole);
167 }
168
169 // the equilibrium enthalpy monotonically increases with T;
170 // if the current value is below the target, then we know the
171 // current temperature is too low. Set the lower bounds.
172 if (Hnow < Htarget) {
173 if (Tnow > Tlow) {
174 Tlow = Tnow;
175 Hlow = Hnow;
176 }
177 } else {
178 // the current enthalpy is greater than the target; therefore
179 // the current temperature is too high. Set the high bounds.
180 if (Tnow < Thigh) {
181 Thigh = Tnow;
182 Hhigh = Hnow;
183 }
184 }
185 double dT;
186 if (Hlow != Undef && Hhigh != Undef) {
187 cpb = (Hhigh - Hlow)/(Thigh - Tlow);
188 dT = (Htarget - Hnow)/cpb;
189 double dTa = fabs(dT);
190 double dTmax = 0.5*fabs(Thigh - Tlow);
191 if (dTa > dTmax) {
192 dT *= dTmax/dTa;
193 }
194 } else {
195 double Tnew = sqrt(Tlow*Thigh);
196 dT = clip(Tnew - Tnow, -200.0, 200.0);
197 }
198 double acpb = std::max(fabs(cpb), 1.0E-6);
199 double denom = std::max(fabs(Htarget), acpb);
200 double Herr = Htarget - Hnow;
201 double HConvErr = fabs((Herr)/denom);
202 if (printLvl > 0) {
203 plogf(" equilibrate_HP: It = %d, Tcurr = %g Hcurr = %g, Htarget = %g\n",
204 n, Tnow, Hnow, Htarget);
205 plogf(" H rel error = %g, cp = %g, HConvErr = %g\n",
206 Herr, cpb, HConvErr);
207 }
208
209 if (HConvErr < err) { // || dTa < 1.0e-4) {
210 if (printLvl > 0) {
211 plogf(" equilibrate_HP: CONVERGENCE: Hfinal = %g Tfinal = %g, Its = %d \n",
212 Hnow, Tnow, n);
213 plogf(" H rel error = %g, cp = %g, HConvErr = %g\n",
214 Herr, cpb, HConvErr);
215 }
216 return iSuccess;
217 }
218 double Tnew = Tnow + dT;
219 if (Tnew < 0.0) {
220 Tnew = 0.5*Tnow;
221 }
222 m_mix->setTemperature(Tnew);
223 } catch (const CanteraError&) {
224 if (!estimateEquil) {
225 strt = -1;
226 } else {
227 double Tnew = 0.5*(Tnow + Thigh);
228 if (fabs(Tnew - Tnow) < 1.0) {
229 Tnew = Tnow + 1.0;
230 }
231 m_mix->setTemperature(Tnew);
232 }
233 }
234 }
235 throw CanteraError("vcs_MultiPhaseEquil::equilibrate_HP",
236 "No convergence for T");
237}
238
239int vcs_MultiPhaseEquil::equilibrate_SP(double Starget, double Tlow, double Thigh,
240 int estimateEquil, int printLvl, double err, int maxsteps, int loglevel)
241{
242 int maxiter = 100;
243 int strt = estimateEquil;
244
245 // Lower bound on T. This will change as we progress in the calculation
246 if (Tlow <= 0.0) {
247 Tlow = 0.5 * m_mix->minTemp();
248 }
249 // Upper bound on T. This will change as we progress in the calculation
250 if (Thigh <= 0.0 || Thigh > 1.0E6) {
251 Thigh = 2.0 * m_mix->maxTemp();
252 }
253
254 double cpb = 1.0, dT;
255 double Slow = Undef;
256 double Shigh = Undef;
257 double Tnow = m_mix->temperature();
258 Tlow = std::min(Tnow, Tlow);
259 Thigh = std::max(Tnow, Thigh);
260 int printLvlSub = std::max(printLvl - 1, 0);
261
262 for (int n = 0; n < maxiter; n++) {
263 // start with a loose error tolerance, but tighten it as we get
264 // close to the final temperature
265 try {
266 Tnow = m_mix->temperature();
267 int iSuccess = equilibrate_TP(strt, printLvlSub, err, maxsteps, loglevel);
268 strt = 0;
269 double Snow = m_mix->entropy();
270 double pmoles[10];
271 pmoles[0] = m_mix->phaseMoles(0);
272 double Tmoles = pmoles[0];
273 double SperMole = Snow/Tmoles;
274 if (printLvl > 0) {
275 plogf("T = %g, Snow = %g ,Tmoles = %g, SperMole = %g\n",
276 Tnow, Snow, Tmoles, SperMole);
277 }
278
279 // the equilibrium entropy monotonically increases with T;
280 // if the current value is below the target, then we know the
281 // current temperature is too low. Set the lower bounds to the
282 // current condition.
283 if (Snow < Starget) {
284 if (Tnow > Tlow) {
285 Tlow = Tnow;
286 Slow = Snow;
287 } else {
288 if (Slow > Starget && Snow < Slow) {
289 Thigh = Tlow;
290 Shigh = Slow;
291 Tlow = Tnow;
292 Slow = Snow;
293 }
294 }
295 } else {
296 // the current enthalpy is greater than the target; therefore
297 // the current temperature is too high. Set the high bounds.
298 if (Tnow < Thigh) {
299 Thigh = Tnow;
300 Shigh = Snow;
301 }
302 }
303 if (Slow != Undef && Shigh != Undef) {
304 cpb = (Shigh - Slow)/(Thigh - Tlow);
305 dT = (Starget - Snow)/cpb;
306 double Tnew = Tnow + dT;
307 double dTa = fabs(dT);
308 double dTmax = 0.5*fabs(Thigh - Tlow);
309 if (Tnew > Thigh || Tnew < Tlow) {
310 dTmax = 1.5*fabs(Thigh - Tlow);
311 }
312 dTmax = std::min(dTmax, 300.);
313 if (dTa > dTmax) {
314 dT *= dTmax/dTa;
315 }
316 } else {
317 double Tnew = sqrt(Tlow*Thigh);
318 dT = Tnew - Tnow;
319 }
320
321 double acpb = std::max(fabs(cpb), 1.0E-6);
322 double denom = std::max(fabs(Starget), acpb);
323 double Serr = Starget - Snow;
324 double SConvErr = fabs((Serr)/denom);
325 if (printLvl > 0) {
326 plogf(" equilibrate_SP: It = %d, Tcurr = %g Scurr = %g, Starget = %g\n",
327 n, Tnow, Snow, Starget);
328 plogf(" S rel error = %g, cp = %g, SConvErr = %g\n",
329 Serr, cpb, SConvErr);
330 }
331
332 if (SConvErr < err) { // || dTa < 1.0e-4) {
333 if (printLvl > 0) {
334 plogf(" equilibrate_SP: CONVERGENCE: Sfinal = %g Tfinal = %g, Its = %d \n",
335 Snow, Tnow, n);
336 plogf(" S rel error = %g, cp = %g, HConvErr = %g\n",
337 Serr, cpb, SConvErr);
338 }
339 return iSuccess;
340 }
341 double Tnew = Tnow + dT;
342 if (Tnew < 0.0) {
343 Tnew = 0.5*Tnow;
344 }
345 m_mix->setTemperature(Tnew);
346 } catch (const CanteraError&) {
347 if (!estimateEquil) {
348 strt = -1;
349 } else {
350 double Tnew = 0.5*(Tnow + Thigh);
351 if (fabs(Tnew - Tnow) < 1.0) {
352 Tnew = Tnow + 1.0;
353 }
354 m_mix->setTemperature(Tnew);
355 }
356 }
357 }
358 throw CanteraError("vcs_MultiPhaseEquil::equilibrate_SP",
359 "No convergence for T");
360}
361
362int vcs_MultiPhaseEquil::equilibrate(int XY, int estimateEquil, int printLvl,
363 double err, int maxsteps, int loglevel)
364{
365 double xtarget;
366 if (XY == TP) {
367 return equilibrate_TP(estimateEquil, printLvl, err, maxsteps, loglevel);
368 } else if (XY == HP || XY == UP) {
369 if (XY == HP) {
370 xtarget = m_mix->enthalpy();
371 } else {
372 xtarget = m_mix->IntEnergy();
373 }
374 double Tlow = 0.5 * m_mix->minTemp();
375 double Thigh = 2.0 * m_mix->maxTemp();
376 return equilibrate_HP(xtarget, XY, Tlow, Thigh,
377 estimateEquil, printLvl, err, maxsteps, loglevel);
378 } else if (XY == SP) {
379 xtarget = m_mix->entropy();
380 double Tlow = 0.5 * m_mix->minTemp();
381 double Thigh = 2.0 * m_mix->maxTemp();
382 return equilibrate_SP(xtarget, Tlow, Thigh,
383 estimateEquil, printLvl, err, maxsteps, loglevel);
384 } else if (XY == TV) {
385 xtarget = m_mix->temperature();
386 return equilibrate_TV(XY, xtarget,
387 estimateEquil, printLvl, err, maxsteps, loglevel);
388 } else if (XY == HV) {
389 xtarget = m_mix->enthalpy();
390 return equilibrate_TV(XY, xtarget,
391 estimateEquil, printLvl, err, maxsteps, loglevel);
392 } else if (XY == UV) {
393 xtarget = m_mix->IntEnergy();
394 return equilibrate_TV(XY, xtarget,
395 estimateEquil, printLvl, err, maxsteps, loglevel);
396 } else if (XY == SV) {
397 xtarget = m_mix->entropy();
398 return equilibrate_TV(XY, xtarget, estimateEquil,
399 printLvl, err, maxsteps, loglevel);
400 } else {
401 throw CanteraError("vcs_MultiPhaseEquil::equilibrate",
402 "Unsupported Option");
403 }
404}
405
406int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil, int printLvl, double err,
407 int maxsteps, int loglevel)
408{
409 int maxit = maxsteps;
410 clockWC tickTock;
411 m_printLvl = printLvl;
412 m_vsolve.m_printLvl = printLvl;
413 m_vsolve.m_doEstimateEquil = estimateEquil;
414
415 // Check obvious bounds on the temperature and pressure NOTE, we may want to
416 // do more here with the real bounds given by the ThermoPhase objects.
417 if (m_mix->temperature() <= 0.0) {
418 throw CanteraError("vcs_MultiPhaseEquil::equilibrate_TP",
419 "Temperature less than zero on input");
420 }
421 if (m_mix->pressure() <= 0.0) {
422 throw CanteraError("vcs_MultiPhaseEquil::equilibrate_TP",
423 "Pressure less than zero on input");
424 }
425
426 //! Call the thermo Program
427 int ip1 = m_printLvl;
428 int ipr = std::max(0, m_printLvl-1);
429 if (m_printLvl >= 3) {
430 ip1 = m_printLvl - 2;
431 } else {
432 ip1 = 0;
433 }
434 int iSuccess = m_vsolve.vcs(ipr, ip1, maxit);
435
436 double te = tickTock.secondsWC();
437 if (printLvl > 0) {
438 vector<double> mu(m_mix->nSpecies());
439 m_mix->getChemPotentials(mu.data());
440 plogf("\n Results from vcs:\n");
441 if (iSuccess != 0) {
442 plogf("\nVCS FAILED TO CONVERGE!\n");
443 }
444 plogf("\n");
445 plogf("Temperature = %g Kelvin\n", m_vsolve.m_temperature);
446 plogf("Pressure = %g Pa\n", m_vsolve.m_pressurePA);
447 plogf("\n");
448 plogf("----------------------------------------"
449 "---------------------\n");
450 plogf(" Name Mole_Number(kmol)");
451 plogf(" Mole_Fraction Chem_Potential (J/kmol)\n");
452 plogf("--------------------------------------------------"
453 "-----------\n");
454 for (size_t i = 0; i < m_mix->nSpecies(); i++) {
455 plogf("%-12s", m_mix->speciesName(i));
457 plogf(" %15.3e %15.3e ", 0.0, m_mix->moleFraction(i));
458 plogf("%15.3e\n", mu[i]);
459 } else {
460 plogf(" %15.3e %15.3e ", m_mix->speciesMoles(i), m_mix->moleFraction(i));
461 if (m_mix->speciesMoles(i) <= 0.0) {
462 size_t iph = m_vsolve.m_phaseID[i];
463 vcs_VolPhase* VPhase = m_vsolve.m_VolPhaseList[iph].get();
464 if (VPhase->nSpecies() > 1) {
465 plogf(" -1.000e+300\n");
466 } else {
467 plogf("%15.3e\n", mu[i]);
468 }
469 } else {
470 plogf("%15.3e\n", mu[i]);
471 }
472 }
473 }
474 plogf("------------------------------------------"
475 "-------------------\n");
476 if (printLvl > 2 && m_vsolve.m_timing_print_lvl > 0) {
477 plogf("Total time = %12.6e seconds\n", te);
478 }
479 }
480 return iSuccess;
481}
482
483void vcs_MultiPhaseEquil::reportCSV(const string& reportFile)
484{
485 size_t nphase = m_vsolve.m_numPhases;
486
487 FILE* FP = fopen(reportFile.c_str(), "w");
488 if (!FP) {
489 throw CanteraError("vcs_MultiPhaseEquil::reportCSV",
490 "Failure to open file");
491 }
492 vector<double> VolPM;
493 vector<double> activity;
494 vector<double> ac;
495 vector<double> mu;
496 vector<double> mu0;
497 vector<double> molalities;
498
499 double vol = 0.0;
500 for (size_t iphase = 0; iphase < nphase; iphase++) {
501 ThermoPhase& tref = m_mix->phase(iphase);
502 size_t nSpecies = tref.nSpecies();
503 VolPM.resize(nSpecies, 0.0);
504 tref.getPartialMolarVolumes(&VolPM[0]);
505 vcs_VolPhase* volP = m_vsolve.m_VolPhaseList[iphase].get();
506
507 double TMolesPhase = volP->totalMoles();
508 double VolPhaseVolumes = 0.0;
509 for (size_t k = 0; k < nSpecies; k++) {
510 VolPhaseVolumes += VolPM[k] * tref.moleFraction(k);
511 }
512 VolPhaseVolumes *= TMolesPhase;
513 vol += VolPhaseVolumes;
514 }
515
516 fprintf(FP,"--------------------- VCS_MULTIPHASE_EQUIL FINAL REPORT"
517 " -----------------------------\n");
518 fprintf(FP,"Temperature = %11.5g kelvin\n", m_mix->temperature());
519 fprintf(FP,"Pressure = %11.5g Pascal\n", m_mix->pressure());
520 fprintf(FP,"Total Volume = %11.5g m**3\n", vol);
521 fprintf(FP,"Number Basis optimizations = %d\n", m_vsolve.m_VCount->Basis_Opts);
522 fprintf(FP,"Number VCS iterations = %d\n", m_vsolve.m_VCount->Its);
523
524 for (size_t iphase = 0; iphase < nphase; iphase++) {
525 ThermoPhase& tref = m_mix->phase(iphase);
526 string phaseName = tref.name();
527 vcs_VolPhase* volP = m_vsolve.m_VolPhaseList[iphase].get();
528 double TMolesPhase = volP->totalMoles();
529 size_t nSpecies = tref.nSpecies();
530 activity.resize(nSpecies, 0.0);
531 ac.resize(nSpecies, 0.0);
532 mu0.resize(nSpecies, 0.0);
533 mu.resize(nSpecies, 0.0);
534 VolPM.resize(nSpecies, 0.0);
535 molalities.resize(nSpecies, 0.0);
536 int actConvention = tref.activityConvention();
537 tref.getActivities(&activity[0]);
538 tref.getActivityCoefficients(&ac[0]);
539 tref.getStandardChemPotentials(&mu0[0]);
540 tref.getPartialMolarVolumes(&VolPM[0]);
541 tref.getChemPotentials(&mu[0]);
542 double VolPhaseVolumes = 0.0;
543 for (size_t k = 0; k < nSpecies; k++) {
544 VolPhaseVolumes += VolPM[k] * tref.moleFraction(k);
545 }
546 VolPhaseVolumes *= TMolesPhase;
547 vol += VolPhaseVolumes;
548
549 if (actConvention == 1) {
550 MolalityVPSSTP* mTP = static_cast<MolalityVPSSTP*>(&tref);
551 mTP->getMolalities(&molalities[0]);
552 tref.getChemPotentials(&mu[0]);
553
554 if (iphase == 0) {
555 fprintf(FP," Name, Phase, PhaseMoles, Mole_Fract, "
556 "Molalities, ActCoeff, Activity,"
557 "ChemPot_SS0, ChemPot, mole_num, PMVol, Phase_Volume\n");
558
559 fprintf(FP," , , (kmol), , "
560 " , , ,"
561 " (J/kmol), (J/kmol), (kmol), (m**3/kmol), (m**3)\n");
562 }
563 for (size_t k = 0; k < nSpecies; k++) {
564 string sName = tref.speciesName(k);
565 fprintf(FP,"%12s, %11s, %11.3e, %11.3e, %11.3e, %11.3e, %11.3e,"
566 "%11.3e, %11.3e, %11.3e, %11.3e, %11.3e\n",
567 sName.c_str(),
568 phaseName.c_str(), TMolesPhase,
569 tref.moleFraction(k), molalities[k], ac[k], activity[k],
570 mu0[k]*1.0E-6, mu[k]*1.0E-6,
571 tref.moleFraction(k) * TMolesPhase,
572 VolPM[k], VolPhaseVolumes);
573 }
574 } else {
575 if (iphase == 0) {
576 fprintf(FP," Name, Phase, PhaseMoles, Mole_Fract, "
577 "Molalities, ActCoeff, Activity,"
578 " ChemPotSS0, ChemPot, mole_num, PMVol, Phase_Volume\n");
579
580 fprintf(FP," , , (kmol), , "
581 " , , ,"
582 " (J/kmol), (J/kmol), (kmol), (m**3/kmol), (m**3)\n");
583 }
584 for (size_t k = 0; k < nSpecies; k++) {
585 molalities[k] = 0.0;
586 }
587 for (size_t k = 0; k < nSpecies; k++) {
588 string sName = tref.speciesName(k);
589 fprintf(FP,"%12s, %11s, %11.3e, %11.3e, %11.3e, %11.3e, %11.3e, "
590 "%11.3e, %11.3e,% 11.3e, %11.3e, %11.3e\n",
591 sName.c_str(),
592 phaseName.c_str(), TMolesPhase,
593 tref.moleFraction(k), molalities[k], ac[k],
594 activity[k], mu0[k]*1.0E-6, mu[k]*1.0E-6,
595 tref.moleFraction(k) * TMolesPhase,
596 VolPM[k], VolPhaseVolumes);
597 }
598 }
599 }
600 fclose(FP);
601}
602
603}
ThermoPhase object for the ideal molal equation of state (see Thermodynamic Properties and class Idea...
Header file for an ideal solid solution model with incompressible thermodynamics (see Thermodynamic P...
Base class for exceptions thrown by Cantera classes.
MolalityVPSSTP is a derived class of ThermoPhase that handles variable pressure standard state method...
void getMolalities(double *const molal) const
This function will return the molalities of the species.
A class for multiphase mixtures.
Definition MultiPhase.h:57
double speciesMoles(size_t kGlob) const
Returns the moles of global species k. units = kmol.
size_t nSpecies() const
Number of species, summed over all phases.
Definition MultiPhase.h:123
double enthalpy() const
The enthalpy of the mixture [J].
double pressure() const
Pressure [Pa].
Definition MultiPhase.h:388
double minTemp() const
Minimum temperature for which all solution phases have valid thermo data.
Definition MultiPhase.h:246
double entropy() const
The entropy of the mixture [J/K].
double temperature() const
Temperature [K].
Definition MultiPhase.h:328
void getChemPotentials(double *mu) const
Returns a vector of Chemical potentials.
double moleFraction(const size_t kGlob) const
Returns the mole fraction of global species k.
void setPressure(double P)
Set the pressure [Pa].
Definition MultiPhase.h:403
string speciesName(const size_t kGlob) const
Name of species with global index kGlob.
double phaseMoles(const size_t n) const
Return the number of moles in phase n.
double volume() const
The total mixture volume [m^3].
void setTemperature(const double T)
Set the temperature [K].
ThermoPhase & phase(size_t n)
Return a reference to phase n.
double maxTemp() const
Maximum temperature for which all solution phases have valid thermo data.
Definition MultiPhase.h:253
double IntEnergy() const
The internal energy of the mixture [J].
size_t nSpecies() const
Returns the number of species in the phase.
Definition Phase.h:231
string speciesName(size_t k) const
Name of the species with index k.
Definition Phase.cpp:142
double moleFraction(size_t k) const
Return the mole fraction of a single species.
Definition Phase.cpp:439
string name() const
Return the name of the phase.
Definition Phase.cpp:20
Base class for a phase with thermodynamic properties.
virtual void getActivityCoefficients(double *ac) const
Get the array of non-dimensional molar-based activity coefficients at the current solution temperatur...
virtual int activityConvention() const
This method returns the convention used in specification of the activities, of which there are curren...
virtual void getActivities(double *a) const
Get the array of non-dimensional activities at the current solution temperature, pressure,...
virtual void getStandardChemPotentials(double *mu) const
Get the array of chemical potentials at unit activity for the species at their standard states at the...
virtual void getChemPotentials(double *mu) const
Get the species chemical potentials. Units: J/kmol.
virtual void getPartialMolarVolumes(double *vbar) const
Return an array of partial molar volumes for the species in the mixture.
int Its
Current number of iterations in the main loop of vcs_TP() to solve for thermo equilibrium.
int Basis_Opts
number of optimizations of the components basis set done
size_t m_numPhases
Number of Phases in the problem.
Definition vcs_solve.h:1069
int m_doEstimateEquil
Setting for whether to do an initial estimate.
Definition vcs_solve.h:1146
vector< size_t > m_phaseID
Mapping from the species number to the phase number.
Definition vcs_solve.h:1358
vector< int > m_speciesUnknownType
Specifies the species unknown type.
Definition vcs_solve.h:1166
int m_timing_print_lvl
printing level of timing information
Definition vcs_solve.h:1506
VCS_COUNTERS * m_VCount
Timing and iteration counters for the vcs object.
Definition vcs_solve.h:1485
int vcs(int ipr, int ip1, int maxit)
Solve an equilibrium problem.
vector< unique_ptr< vcs_VolPhase > > m_VolPhaseList
Array of Phase Structures. Length = number of phases.
Definition vcs_solve.h:1395
double m_pressurePA
Pressure.
Definition vcs_solve.h:1274
int m_printLvl
Print level for print routines.
Definition vcs_solve.h:985
double m_temperature
Temperature (Kelvin)
Definition vcs_solve.h:1271
The class provides the wall clock timer in seconds.
Definition clockWC.h:47
double secondsWC()
Returns the wall clock time in seconds since the last reset.
Definition clockWC.cpp:32
VCS_SOLVE m_vsolve
The object that contains the problem statement and does all of the equilibration work.
int equilibrate_SP(double Starget, double Tlow, double Thigh, int estimateEquil=0, int printLvl=0, double err=1.0E-6, int maxsteps=VCS_MAXSTEPS, int loglevel=-99)
Equilibrate the solution using the current element abundances stored in the MultiPhase object using c...
void reportCSV(const string &reportFile)
Report the equilibrium answer in a comma separated table format.
vcs_MultiPhaseEquil(MultiPhase *mix, int printLvl)
Constructor for the multiphase equilibrium solver.
int equilibrate_TP(int estimateEquil=0, int printLvl=0, double err=1.0e-6, int maxsteps=VCS_MAXSTEPS, int loglevel=-99)
Equilibrate the solution using the current element abundances stored in the MultiPhase object using c...
int equilibrate(int XY, int estimateEquil=0, int printLvl=0, double err=1.0e-6, int maxsteps=VCS_MAXSTEPS, int loglevel=-99)
Equilibrate the solution using the current element abundances stored in the MultiPhase object.
int equilibrate_TV(int XY, double xtarget, int estimateEquil=0, int printLvl=0, double err=1.0E-6, int maxsteps=VCS_MAXSTEPS, int logLevel=-99)
Equilibrate the solution using the current element abundances stored in the MultiPhase object using c...
int m_printLvl
Print level from the VCSnonlinear package.
int equilibrate_HP(double Htarget, int XY, double Tlow, double Thigh, int estimateEquil=0, int printLvl=0, double err=1.0E-6, int maxsteps=VCS_MAXSTEPS, int loglevel=-99)
Equilibrate the solution using the current element abundances stored in the MultiPhase object using e...
MultiPhase * m_mix
Pointer to the MultiPhase mixture that will be equilibrated.
Phase information and Phase calculations for vcs.
size_t nSpecies() const
Return the number of species in the phase.
double totalMoles() const
Return the total moles in the phase.
Declarations for a simple class that implements an Ansi C wall clock timer (see Cantera::clockWC).
T clip(const T &value, const T &lower, const T &upper)
Clip value such that lower <= value <= upper.
Definition global.h:329
Namespace for the Cantera kernel.
Definition AnyMap.cpp:564
const double Undef
Fairly random number to be used to initialize variables against to see if they are subsequently defin...
Definition ct_defs.h:164
Contains const definitions for types of species reference-state thermodynamics managers (see Species ...
Contains declarations for string manipulation functions within Cantera.
Interface class for the vcsnonlinear solver.
Header for the object representing each phase within vcs.
#define VCS_SPECIES_TYPE_INTERFACIALVOLTAGE
Unknown refers to the voltage level of a phase.
Definition vcs_defs.h:286
#define plogf
define this Cantera function to replace printf