Cantera  3.1.0a1
Loading...
Searching...
No Matches
PDSS_HKFT.cpp
Go to the documentation of this file.
1/**
2 * @file PDSS_HKFT.cpp
3 * Definitions for the class PDSS_HKFT (pressure dependent standard state)
4 * which handles calculations for a single species in a phase using the
5 * HKFT standard state
6 * (see @ref pdssthermo and class @link Cantera::PDSS_HKFT PDSS_HKFT@endlink).
7 */
8
9// This file is part of Cantera. See License.txt in the top-level directory or
10// at https://cantera.org/license.txt for license and copyright information.
11
16#include "cantera/base/global.h"
17
18namespace Cantera
19{
20// Set the default to error exit if there is an input file inconsistency
22
24{
25 m_pres = OneAtm;
26 m_units.setDefaults({"cm", "gmol", "cal", "bar"});
27}
28
30{
31 // Ok we may change this evaluation method in the future.
32 double h = gibbs_mole() + m_temp * entropy_mole();
33 return h;
34}
35
37{
38 return enthalpy_RT() - molarVolume() * m_pres;
39}
40
42{
43 return m_units.convertTo(m_Entrop_tr_pr, "J/kmol") + deltaS();
44}
45
47{
48 return m_Mu0_tr_pr + deltaG();
49}
50
51double PDSS_HKFT::cp_mole() const
52{
53 double pbar = m_pres * 1.0E-5;
54 double c1term = m_c1;
55 double c2term = m_c2 / (m_temp - 228.) / (m_temp - 228.);
56 double a3term = -m_a3 / (m_temp - 228.) / (m_temp - 228.) / (m_temp - 228.) * 2.0 * m_temp * (pbar - m_presR_bar);
57 double a4term = -m_a4 / (m_temp - 228.) / (m_temp - 228.) / (m_temp - 228.) * 2.0 * m_temp
58 * log((2600. + pbar)/(2600. + m_presR_bar));
59
60 double omega_j;
61 double domega_jdT;
62 double d2omega_jdT2;
63 if (m_charge_j == 0.0) {
64 omega_j = m_omega_pr_tr;
65 domega_jdT = 0.0;
66 d2omega_jdT2 = 0.0;
67 } else {
68 double nu = 166027;
69 double r_e_j_pr_tr = m_charge_j * m_charge_j / (m_omega_pr_tr/nu + m_charge_j/3.082);
70 double gval = gstar(m_temp, m_pres, 0);
71 double dgvaldT = gstar(m_temp, m_pres, 1);
72 double d2gvaldT2 = gstar(m_temp, m_pres, 2);
73
74 double r_e_j = r_e_j_pr_tr + fabs(m_charge_j) * gval;
75 double dr_e_jdT = fabs(m_charge_j) * dgvaldT;
76 double d2r_e_jdT2 = fabs(m_charge_j) * d2gvaldT2;
77 double r_e_j2 = r_e_j * r_e_j;
78
79 double charge2 = m_charge_j * m_charge_j;
80 double r_e_H = 3.082 + gval;
81 double r_e_H2 = r_e_H * r_e_H;
82 omega_j = nu * (charge2 / r_e_j - m_charge_j / r_e_H);
83 domega_jdT = nu * (-(charge2 / r_e_j2 * dr_e_jdT)
84 +(m_charge_j / r_e_H2 * dgvaldT));
85 d2omega_jdT2 = nu * (2.0*charge2*dr_e_jdT*dr_e_jdT/(r_e_j2*r_e_j) - charge2*d2r_e_jdT2/r_e_j2
86 -2.0*m_charge_j*dgvaldT*dgvaldT/(r_e_H2*r_e_H) + m_charge_j*d2gvaldT2 /r_e_H2);
87 }
88
89 double relepsilon = m_waterProps->relEpsilon(m_temp, m_pres, 0);
90 double drelepsilondT = m_waterProps->relEpsilon(m_temp, m_pres, 1);
91 double Y = drelepsilondT / (relepsilon * relepsilon);
92 double d2relepsilondT2 = m_waterProps->relEpsilon(m_temp, m_pres, 2);
93
94 double X = d2relepsilondT2 / (relepsilon* relepsilon) - 2.0 * relepsilon * Y * Y;
95 double Z = -1.0 / relepsilon;
96
97 double yterm = 2.0 * m_temp * Y * domega_jdT;
98 double xterm = omega_j * m_temp * X;
99 double otterm = m_temp * d2omega_jdT2 * (Z + 1.0);
100 double rterm = - m_domega_jdT_prtr * (m_Z_pr_tr + 1.0);
101
102 double Cp_calgmol = c1term + c2term + a3term + a4term + yterm + xterm + otterm + rterm;
103
104 return m_units.convertTo(Cp_calgmol, "J/kmol");
105}
106
108{
109 // Initially do all calculations in (cal/gmol/Pa)
110
111 double a1term = m_a1 * 1.0E-5;
112 double a2term = m_a2 / (2600.E5 + m_pres);
113 double a3term = m_a3 * 1.0E-5/ (m_temp - 228.);
114 double a4term = m_a4 / (m_temp - 228.) / (2600.E5 + m_pres);
115
116 double omega_j;
117 double domega_jdP;
118 if (m_charge_j == 0.0) {
119 omega_j = m_omega_pr_tr;
120 domega_jdP = 0.0;
121 } else {
122 double nu = 166027.;
123 double charge2 = m_charge_j * m_charge_j;
124 double r_e_j_pr_tr = charge2 / (m_omega_pr_tr/nu + m_charge_j/3.082);
125
126 double gval = gstar(m_temp, m_pres, 0);
127 double dgvaldP = gstar(m_temp, m_pres, 3);
128
129 double r_e_j = r_e_j_pr_tr + fabs(m_charge_j) * gval;
130 double r_e_H = 3.082 + gval;
131
132 omega_j = nu * (charge2 / r_e_j - m_charge_j / r_e_H);
133 double dr_e_jdP = fabs(m_charge_j) * dgvaldP;
134 domega_jdP = - nu * (charge2 / (r_e_j * r_e_j) * dr_e_jdP)
135 + nu * m_charge_j / (r_e_H * r_e_H) * dgvaldP;
136 }
137
138 double drelepsilondP = m_waterProps->relEpsilon(m_temp, m_pres, 3);
139 double relepsilon = m_waterProps->relEpsilon(m_temp, m_pres, 0);
140 double Q = drelepsilondP / (relepsilon * relepsilon);
141 double Z = -1.0 / relepsilon;
142 double wterm = - domega_jdP * (Z + 1.0);
143 double qterm = - omega_j * Q;
144 double molVol_calgmolPascal = a1term + a2term + a3term + a4term + wterm + qterm;
145
146 // Convert to m**3 / kmol from (cal/gmol/Pa)
147 return m_units.convertTo(molVol_calgmolPascal, "J/kmol");
148}
149
150double PDSS_HKFT::density() const
151{
152 return m_mw / molarVolume();
153}
154
156{
157 double m_psave = m_pres;
159 double ee = gibbs_RT();
160 m_pres = m_psave;
161 return ee;
162}
163
165{
166 double m_psave = m_pres;
168 double hh = enthalpy_RT();
169 m_pres = m_psave;
170 return hh;
171}
172
174{
175 double m_psave = m_pres;
177 double ee = entropy_R();
178 m_pres = m_psave;
179 return ee;
180}
181
183{
184 double m_psave = m_pres;
186 double ee = cp_R();
187 m_pres = m_psave;
188 return ee;
189}
190
192{
193 double m_psave = m_pres;
195 double ee = molarVolume();
196 m_pres = m_psave;
197 return ee;
198}
199
200void PDSS_HKFT::setState_TP(double temp, double pres)
201{
202 setTemperature(temp);
203 setPressure(pres);
204}
205
207{
209 if (m_input.hasKey("h0")) {
210 m_deltaH_formation_tr_pr = m_input.convert("h0", "cal/gmol");
211 }
212 if (m_input.hasKey("g0")) {
213 m_deltaG_formation_tr_pr = m_input.convert("g0", "cal/gmol");
214 }
215 if (m_input.hasKey("s0")) {
216 m_Entrop_tr_pr = m_input.convert("s0", "cal/gmol/K");
217 }
218 auto& units = m_input.units();
219 if (m_input.hasKey("a")) {
220 auto& a = m_input["a"].asVector<AnyValue>(4);
221 m_a1 = units.convert(a[0], "cal/gmol/bar");
222 m_a2 = units.convert(a[1], "cal/gmol");
223 m_a3 = units.convert(a[2], "cal*K/gmol/bar");
224 m_a4 = units.convert(a[3], "cal*K/gmol");
225 }
226 if (m_input.hasKey("c")) {
227 auto& c = m_input["c"].asVector<AnyValue>(2);
228 m_c1 = units.convert(c[0], "cal/gmol/K");
229 m_c2 = units.convert(c[1], "cal*K/gmol");
230 }
231 if (m_input.hasKey("omega")) {
232 m_omega_pr_tr = m_input.convert("omega", "cal/gmol");
233 }
234
235 // Ok, if we are missing one, then we construct its value from the other two.
236 // This code has been internally verified.
238 if (std::isnan(m_deltaH_formation_tr_pr)) {
240 double Hcalc = m_Mu0_tr_pr + 298.15 * m_units.convertTo(m_Entrop_tr_pr, "J/kmol");
241 m_deltaH_formation_tr_pr = m_units.convertFrom(Hcalc, "J/kmol");
242 } else if (std::isnan(m_deltaG_formation_tr_pr)) {
243 double DHjmol = m_units.convertTo(m_deltaH_formation_tr_pr, "J/kmol");
244 m_Mu0_tr_pr = DHjmol - 298.15 * m_units.convertTo(m_Entrop_tr_pr, "J/kmol");
246 double tmp = m_Mu0_tr_pr;
248 double totalSum = m_Mu0_tr_pr - tmp;
249 m_Mu0_tr_pr = tmp;
250 m_deltaG_formation_tr_pr = m_units.convertFrom(m_Mu0_tr_pr - totalSum, "J/kmol");
251 } else if (std::isnan(m_Entrop_tr_pr)) {
253 double DHjmol = m_units.convertTo(m_deltaH_formation_tr_pr, "J/kmol");
254 m_Entrop_tr_pr = m_units.convertFrom((DHjmol - m_Mu0_tr_pr) / 298.15, "J/kmol");
255 }
256
257 m_waterSS = &dynamic_cast<PDSS_Water&>(*m_tp->providePDSS(0));
258
259 // Section to initialize m_Z_pr_tr and m_Y_pr_tr
260 m_temp = 273.15 + 25.;
261 m_pres = OneAtm;
262 double relepsilon = m_waterProps->relEpsilon(m_temp, m_pres, 0);
265 m_Z_pr_tr = -1.0 / relepsilon;
266 double drelepsilondT = m_waterProps->relEpsilon(m_temp, m_pres, 1);
267 m_Y_pr_tr = drelepsilondT / (relepsilon * relepsilon);
268 m_waterProps = make_unique<WaterProps>(m_waterSS);
269 m_presR_bar = OneAtm / 1.0E5;
270 m_presR_bar = 1.0;
272
273 // Ok, we have mu. Let's check it against the input value
274 // of DH_F to see that we have some internal consistency
275 double Hcalc = m_Mu0_tr_pr + 298.15 * m_units.convertTo(m_Entrop_tr_pr, "J/kmol");
276 double DHjmol = m_units.convertTo(m_deltaH_formation_tr_pr, "J/kmol");
277
278 // If the discrepancy is greater than 100 cal gmol-1, print
279 // an error and exit.
280 if (fabs(Hcalc -DHjmol) > m_units.convertTo(100, "J/kmol")) {
281 string sname = m_tp->speciesName(m_spindex);
283 throw CanteraError("PDSS_HKFT::initThermo", "For {}, DHjmol is"
284 " not consistent with G and S: {} vs {} J/kmol",
285 sname, Hcalc, m_units.convertTo(m_deltaH_formation_tr_pr, "J/kmol"));
286 } else {
287 warn_user("PDSS_HKFT::initThermo",
288 "DHjmol for {} is not consistent with G and S: calculated {} "
289 "vs input {} J/kmol; continuing with consistent DHjmol = {}",
290 sname, Hcalc, m_units.convertTo(m_deltaH_formation_tr_pr, "J/kmol"),
291 Hcalc);
292 m_deltaH_formation_tr_pr = m_units.convertFrom(Hcalc, "J/kmol");
293 }
294 }
295 double nu = 166027;
296 double r_e_j_pr_tr;
297 if (m_charge_j != 0.0) {
298 r_e_j_pr_tr = m_charge_j * m_charge_j / (m_omega_pr_tr/nu + m_charge_j/3.082);
299 } else {
300 r_e_j_pr_tr = 0.0;
301 }
302
303 if (m_charge_j == 0.0) {
304 m_domega_jdT_prtr = 0.0;
305 } else {
306 double gval = gstar(m_temp, m_pres, 0);
307 double dgvaldT = gstar(m_temp, m_pres, 1);
308 double r_e_j = r_e_j_pr_tr + fabs(m_charge_j) * gval;
309 double dr_e_jdT = fabs(m_charge_j) * dgvaldT;
310 m_domega_jdT_prtr = - nu * (m_charge_j * m_charge_j / (r_e_j * r_e_j) * dr_e_jdT)
311 + nu * m_charge_j / (3.082 + gval) / (3.082 + gval) * dgvaldT;
312 }
313}
314
315void PDSS_HKFT::setDeltaH0(double dh0) {
316 m_deltaH_formation_tr_pr = m_units.convertFrom(dh0, "J/kmol");
317}
318
319void PDSS_HKFT::setDeltaG0(double dg0) {
320 m_deltaG_formation_tr_pr = m_units.convertFrom(dg0, "J/kmol");
321}
322
323void PDSS_HKFT::setS0(double s0) {
324 m_Entrop_tr_pr = m_units.convertFrom(s0, "J/kmol/K");
325}
326
327void PDSS_HKFT::set_a(double* a) {
328 m_a1 = m_units.convertFrom(a[0], "J/kmol/Pa");
329 m_a2 = m_units.convertFrom(a[1], "J/kmol");
330 m_a3 = m_units.convertFrom(a[2], "J*K/kmol/Pa");
331 m_a4 = m_units.convertFrom(a[3], "J*K/kmol");
332}
333
334void PDSS_HKFT::set_c(double* c) {
335 m_c1 = m_units.convertFrom(c[0], "J/kmol/K");
336 m_c2 = m_units.convertFrom(c[1], "J*K/kmol");
337}
338
339void PDSS_HKFT::setOmega(double omega) {
340 m_omega_pr_tr = m_units.convertFrom(omega, "J/kmol");
341}
342
344{
345 PDSS::getParameters(eosNode);
346 eosNode["model"] = "HKFT";
347 eosNode["h0"].setQuantity(m_deltaH_formation_tr_pr, "cal/gmol");
348 eosNode["g0"].setQuantity(m_deltaG_formation_tr_pr, "cal/gmol");
349 eosNode["s0"].setQuantity(m_Entrop_tr_pr, "cal/gmol/K");
350
351 vector<AnyValue> a(4), c(2);
352 a[0].setQuantity(m_a1, "cal/gmol/bar");
353 a[1].setQuantity(m_a2, "cal/gmol");
354 a[2].setQuantity(m_a3, "cal*K/gmol/bar");
355 a[3].setQuantity(m_a4, "cal*K/gmol");
356 eosNode["a"] = std::move(a);
357
358 c[0].setQuantity(m_c1, "cal/gmol/K");
359 c[1].setQuantity(m_c2, "cal*K/gmol");
360 eosNode["c"] = std::move(c);
361
362 eosNode["omega"].setQuantity(m_omega_pr_tr, "cal/gmol");
363}
364
365double PDSS_HKFT::deltaG() const
366{
367 double pbar = m_pres * 1.0E-5;
368 double sterm = - m_Entrop_tr_pr * (m_temp - 298.15);
369 double c1term = -m_c1 * (m_temp * log(m_temp/298.15) - (m_temp - 298.15));
370 double a1term = m_a1 * (pbar - m_presR_bar);
371 double a2term = m_a2 * log((2600. + pbar)/(2600. + m_presR_bar));
372 double c2term = -m_c2 * ((1.0/(m_temp - 228.) - 1.0/(298.15 - 228.)) * (228. - m_temp)/228.
373 - m_temp / (228.*228.) * log((298.15*(m_temp-228.)) / (m_temp*(298.15-228.))));
374 double a3term = m_a3 / (m_temp - 228.) * (pbar - m_presR_bar);
375 double a4term = m_a4 / (m_temp - 228.) * log((2600. + pbar)/(2600. + m_presR_bar));
376
377 double omega_j;
378 if (m_charge_j == 0.0) {
379 omega_j = m_omega_pr_tr;
380 } else {
381 double nu = 166027;
382 double r_e_j_pr_tr = m_charge_j * m_charge_j / (m_omega_pr_tr/nu + m_charge_j/3.082);
383 double gval = gstar(m_temp, m_pres, 0);
384 double r_e_j = r_e_j_pr_tr + fabs(m_charge_j) * gval;
385 omega_j = nu * (m_charge_j * m_charge_j / r_e_j - m_charge_j / (3.082 + gval));
386 }
387
388 double relepsilon = m_waterProps->relEpsilon(m_temp, m_pres, 0);
389 double Z = -1.0 / relepsilon;
390 double wterm = - omega_j * (Z + 1.0);
391 double wrterm = m_omega_pr_tr * (m_Z_pr_tr + 1.0);
392 double yterm = m_omega_pr_tr * m_Y_pr_tr * (m_temp - 298.15);
393 double deltaG_calgmol = sterm + c1term + a1term + a2term + c2term + a3term + a4term + wterm + wrterm + yterm;
394
395 return m_units.convertTo(deltaG_calgmol, "J/kmol");
396}
397
398double PDSS_HKFT::deltaS() const
399{
400 double pbar = m_pres * 1.0E-5;
401
402 double c1term = m_c1 * log(m_temp/298.15);
403 double c2term = -m_c2 / 228. * ((1.0/(m_temp - 228.) - 1.0/(298.15 - 228.))
404 + 1.0 / 228. * log((298.15*(m_temp-228.)) / (m_temp*(298.15-228.))));
405 double a3term = m_a3 / (m_temp - 228.) / (m_temp - 228.) * (pbar - m_presR_bar);
406 double a4term = m_a4 / (m_temp - 228.) / (m_temp - 228.) * log((2600. + pbar)/(2600. + m_presR_bar));
407
408 double omega_j;
409 double domega_jdT;
410 if (m_charge_j == 0.0) {
411 omega_j = m_omega_pr_tr;
412 domega_jdT = 0.0;
413 } else {
414 double nu = 166027;
415 double r_e_j_pr_tr = m_charge_j * m_charge_j / (m_omega_pr_tr/nu + m_charge_j/3.082);
416 double gval = gstar(m_temp, m_pres, 0);
417 double dgvaldT = gstar(m_temp, m_pres, 1);
418 double r_e_j = r_e_j_pr_tr + fabs(m_charge_j) * gval;
419 double dr_e_jdT = fabs(m_charge_j) * dgvaldT;
420 omega_j = nu * (m_charge_j * m_charge_j / r_e_j - m_charge_j / (3.082 + gval));
421 domega_jdT = - nu * (m_charge_j * m_charge_j / (r_e_j * r_e_j) * dr_e_jdT)
422 + nu * m_charge_j / (3.082 + gval) / (3.082 + gval) * dgvaldT;
423 }
424
425 double relepsilon = m_waterProps->relEpsilon(m_temp, m_pres, 0);
426 double drelepsilondT = m_waterProps->relEpsilon(m_temp, m_pres, 1);
427 double Y = drelepsilondT / (relepsilon * relepsilon);
428 double Z = -1.0 / relepsilon;
429 double wterm = omega_j * Y;
430 double wrterm = - m_omega_pr_tr * m_Y_pr_tr;
431 double otterm = domega_jdT * (Z + 1.0);
432 double otrterm = - m_domega_jdT_prtr * (m_Z_pr_tr + 1.0);
433 double deltaS_calgmol = c1term + c2term + a3term + a4term + wterm + wrterm + otterm + otrterm;
434
435 return m_units.convertTo(deltaS_calgmol, "J/kmol/K");
436}
437
438double PDSS_HKFT::ag(const double temp, const int ifunc) const
439{
440 static double ag_coeff[3] = { -2.037662, 5.747000E-3, -6.557892E-6};
441 if (ifunc == 0) {
442 return ag_coeff[0] + ag_coeff[1] * temp + ag_coeff[2] * temp * temp;
443 } else if (ifunc == 1) {
444 return ag_coeff[1] + ag_coeff[2] * 2.0 * temp;
445 }
446 if (ifunc != 2) {
447 return 0.0;
448 }
449 return ag_coeff[2] * 2.0;
450}
451
452double PDSS_HKFT::bg(const double temp, const int ifunc) const
453{
454 static double bg_coeff[3] = { 6.107361, -1.074377E-2, 1.268348E-5};
455 if (ifunc == 0) {
456 return bg_coeff[0] + bg_coeff[1] * temp + bg_coeff[2] * temp * temp;
457 } else if (ifunc == 1) {
458 return bg_coeff[1] + bg_coeff[2] * 2.0 * temp;
459 }
460 if (ifunc != 2) {
461 return 0.0;
462 }
463 return bg_coeff[2] * 2.0;
464}
465
466double PDSS_HKFT::f(const double temp, const double pres, const int ifunc) const
467{
468 static double af_coeff[3] = { 3.666666E1, -0.1504956E-9, 0.5107997E-13};
469 double TC = temp - 273.15;
470 double presBar = pres / 1.0E5;
471
472 if (TC < 155.0) {
473 return 0.0;
474 }
475 TC = std::min(TC, 355.0);
476 if (presBar > 1000.) {
477 return 0.0;
478 }
479
480 double T1 = (TC-155.0)/300.;
481 double p2 = (1000. - presBar) * (1000. - presBar);
482 double p3 = (1000. - presBar) * p2;
483 double p4 = p2 * p2;
484 double fac2 = af_coeff[1] * p3 + af_coeff[2] * p4;
485 if (ifunc == 0) {
486 return pow(T1,4.8) + af_coeff[0] * pow(T1, 16.0) * fac2;
487 } else if (ifunc == 1) {
488 return (4.8 * pow(T1,3.8) + 16.0 * af_coeff[0] * pow(T1, 15.0)) / 300. * fac2;
489 } else if (ifunc == 2) {
490 return (4.8 * 3.8 * pow(T1,2.8) + 16.0 * 15.0 * af_coeff[0] * pow(T1, 14.0)) / (300. * 300.) * fac2;
491 } else if (ifunc == 3) {
492 double fac1 = pow(T1,4.8) + af_coeff[0] * pow(T1, 16.0);
493 fac2 = - (3.0 * af_coeff[1] * p2 + 4.0 * af_coeff[2] * p3)/ 1.0E5;
494 return fac1 * fac2;
495 } else {
496 throw NotImplementedError("PDSS_HKFT::f");
497 }
498}
499
500double PDSS_HKFT::g(const double temp, const double pres, const int ifunc) const
501{
502 double afunc = ag(temp, 0);
503 double bfunc = bg(temp, 0);
504 m_waterSS->setState_TP(temp, pres);
506 // density in gm cm-3
507 double dens = m_densWaterSS * 1.0E-3;
508 double gval = afunc * pow((1.0-dens), bfunc);
509 if (dens >= 1.0) {
510 return 0.0;
511 }
512 if (ifunc == 0) {
513 return gval;
514 } else if (ifunc == 1 || ifunc == 2) {
515 double afuncdT = ag(temp, 1);
516 double bfuncdT = bg(temp, 1);
517 double alpha = m_waterSS->thermalExpansionCoeff();
518
519 double fac1 = afuncdT * gval / afunc;
520 double fac2 = bfuncdT * gval * log(1.0 - dens);
521 double fac3 = gval * alpha * bfunc * dens / (1.0 - dens);
522
523 double dgdt = fac1 + fac2 + fac3;
524 if (ifunc == 1) {
525 return dgdt;
526 }
527
528 double afuncdT2 = ag(temp, 2);
529 double bfuncdT2 = bg(temp, 2);
530 double dfac1dT = dgdt * afuncdT / afunc + afuncdT2 * gval / afunc
531 - afuncdT * afuncdT * gval / (afunc * afunc);
532 double ddensdT = - alpha * dens;
533 double dfac2dT = bfuncdT2 * gval * log(1.0 - dens)
534 + bfuncdT * dgdt * log(1.0 - dens)
535 - bfuncdT * gval /(1.0 - dens) * ddensdT;
536 double dalphadT = m_waterSS->dthermalExpansionCoeffdT();
537 double dfac3dT = dgdt * alpha * bfunc * dens / (1.0 - dens)
538 + gval * dalphadT * bfunc * dens / (1.0 - dens)
539 + gval * alpha * bfuncdT * dens / (1.0 - dens)
540 + gval * alpha * bfunc * ddensdT / (1.0 - dens)
541 + gval * alpha * bfunc * dens / ((1.0 - dens) * (1.0 - dens)) * ddensdT;
542
543 return dfac1dT + dfac2dT + dfac3dT;
544 } else if (ifunc == 3) {
545 double beta = m_waterSS->isothermalCompressibility();
546 return - bfunc * gval * dens * beta / (1.0 - dens);
547 } else {
548 throw NotImplementedError("PDSS_HKFT::g");
549 }
550 return 0.0;
551}
552
553double PDSS_HKFT::gstar(const double temp, const double pres, const int ifunc) const
554{
555 double gval = g(temp, pres, ifunc);
556 double fval = f(temp, pres, ifunc);
557 double res = gval - fval;
558 return res;
559}
560
561double PDSS_HKFT::LookupGe(const string& elemName)
562{
563 size_t iE = m_tp->elementIndex(elemName);
564 if (iE == npos) {
565 throw CanteraError("PDSS_HKFT::LookupGe", "element " + elemName + " not found");
566 }
567 double geValue = m_tp->entropyElement298(iE);
568 if (geValue == ENTROPY298_UNKNOWN) {
569 throw CanteraError("PDSS_HKFT::LookupGe",
570 "element " + elemName + " does not have a supplied entropy298");
571 }
572 return geValue * -298.15;
573}
574
576{
577 // Ok let's get the element compositions and conversion factors.
578 double totalSum = 0.0;
579 for (size_t m = 0; m < m_tp->nElements(); m++) {
580 double na = m_tp->nAtoms(m_spindex, m);
581 if (na > 0.0) {
582 totalSum += na * LookupGe(m_tp->elementName(m));
583 }
584 }
585 // Add in the charge
586 if (m_charge_j != 0.0) {
587 totalSum -= m_charge_j * LookupGe("H");
588 }
589 double dg = m_units.convertTo(m_deltaG_formation_tr_pr, "J/kmol");
590 //! Store the result into an internal variable.
591 m_Mu0_tr_pr = dg + totalSum;
592}
593
594}
#define ENTROPY298_UNKNOWN
Number indicating we don't know the entropy of the element in its most stable state at 298....
Definition Elements.h:85
Declarations for the class PDSS_HKFT (pressure dependent standard state) which handles calculations f...
Implementation of a pressure dependent standard state virtual function for a Pure Water Phase (see Sp...
Header file for a derived class of ThermoPhase that handles variable pressure standard state methods ...
A map of string keys to values whose type can vary at runtime.
Definition AnyMap.h:427
bool hasKey(const string &key) const
Returns true if the map contains an item named key.
Definition AnyMap.cpp:1423
const UnitSystem & units() const
Return the default units that should be used to convert stored values.
Definition AnyMap.h:630
double convert(const string &key, const string &units) const
Convert the item stored by the given key to the units specified in units.
Definition AnyMap.cpp:1535
A wrapper for a variable whose type is determined at runtime.
Definition AnyMap.h:86
Base class for exceptions thrown by Cantera classes.
An error indicating that an unimplemented function has been called.
void setOmega(double omega)
Set omega [J/kmol].
double molarVolume() const override
Return the molar volume at standard state.
static int s_InputInconsistencyErrorExit
Static variable determining error exiting.
Definition PDSS_HKFT.h:297
double m_Y_pr_tr
y = dZdT = 1/(esp*esp) desp/dT at 298.15 and 1 bar
Definition PDSS_HKFT.h:278
double m_a4
Input a4 coefficient (cal K gmol-1)
Definition PDSS_HKFT.h:266
double enthalpy_mole() const override
Return the molar enthalpy in units of J kmol-1.
Definition PDSS_HKFT.cpp:29
unique_ptr< WaterProps > m_waterProps
Pointer to the water property calculator.
Definition PDSS_HKFT.h:221
double deltaS() const
Main routine that actually calculates the entropy difference between the reference state at Tr,...
double f(const double temp, const double pres, const int ifunc=0) const
Difference function f appearing in the formulation.
size_t m_spindex
Index of this species within the parent phase.
Definition PDSS_HKFT.h:101
double m_densWaterSS
density of standard-state water. internal temporary variable
Definition PDSS_HKFT.h:218
double LookupGe(const string &elemName)
Function to look up Element Free Energies.
double gibbs_RT_ref() const override
Return the molar Gibbs free energy divided by RT at reference pressure.
VPStandardStateTP * m_tp
Parent VPStandardStateTP (ThermoPhase) object.
Definition PDSS_HKFT.h:100
double m_Entrop_tr_pr
Input value of S_j at Tr and Pr (cal gmol-1 K-1)
Definition PDSS_HKFT.h:254
double m_deltaG_formation_tr_pr
Input value of deltaG of Formation at Tr and Pr (cal gmol-1)
Definition PDSS_HKFT.h:230
void setDeltaG0(double dg0)
Set Gibbs free energy of formation at Pr, Tr [J/kmol].
void initThermo() override
Initialization routine.
double m_a3
Input a3 coefficient (cal K gmol-1 bar-1)
Definition PDSS_HKFT.h:263
double gstar(const double temp, const double pres, const int ifunc=0) const
Evaluate the Gstar value appearing in the HKFT formulation.
double m_charge_j
Charge of the ion.
Definition PDSS_HKFT.h:290
PDSS_HKFT()
Default Constructor.
Definition PDSS_HKFT.cpp:23
double entropy_R_ref() const override
Return the molar entropy divided by R at reference pressure.
double deltaG() const
Main routine that actually calculates the Gibbs free energy difference between the reference state at...
double enthalpy_RT_ref() const override
Return the molar enthalpy divided by RT at reference pressure.
void setS0(double s0)
Set entropy of formation at Pr, Tr [J/kmol/K].
double bg(const double temp, const int ifunc=0) const
Internal formula for the calculation of b_g()
void getParameters(AnyMap &eosNode) const override
Store the parameters needed to reconstruct a copy of this PDSS object.
double m_c2
Input c2 coefficient (cal K gmol-1)
Definition PDSS_HKFT.h:272
double m_presR_bar
Reference pressure is 1 atm in units of bar= 1.0132.
Definition PDSS_HKFT.h:284
void set_a(double *a)
Set "a" coefficients (array of 4 elements).
double intEnergy_mole() const override
Return the molar internal Energy in units of J kmol-1.
Definition PDSS_HKFT.cpp:36
double m_domega_jdT_prtr
small value that is not quite zero
Definition PDSS_HKFT.h:287
double entropy_mole() const override
Return the molar entropy in units of J kmol-1 K-1.
Definition PDSS_HKFT.cpp:41
double m_omega_pr_tr
Input omega_pr_tr coefficient(cal gmol-1)
Definition PDSS_HKFT.h:275
double g(const double temp, const double pres, const int ifunc=0) const
function g appearing in the formulation
void setState_TP(double temp, double pres) override
Set the internal temperature and pressure.
double m_Z_pr_tr
Z = -1 / relEpsilon at 298.15 and 1 bar.
Definition PDSS_HKFT.h:281
void convertDGFormation()
Translate a Gibbs free energy of formation value to a NIST-based Chemical potential.
double cp_mole() const override
Return the molar const pressure heat capacity in units of J kmol-1 K-1.
Definition PDSS_HKFT.cpp:51
double m_deltaH_formation_tr_pr
Input value of deltaH of Formation at Tr and Pr (cal gmol-1)
Definition PDSS_HKFT.h:239
double m_a1
Input a1 coefficient (cal gmol-1 bar-1)
Definition PDSS_HKFT.h:257
double ag(const double temp, const int ifunc=0) const
Internal formula for the calculation of a_g()
double density() const override
Return the standard state density at standard state.
double gibbs_mole() const override
Return the molar Gibbs free energy in units of J kmol-1.
Definition PDSS_HKFT.cpp:46
double m_c1
Input c1 coefficient (cal gmol-1 K-1)
Definition PDSS_HKFT.h:269
double molarVolume_ref() const override
Return the molar volume at reference pressure.
PDSS_Water * m_waterSS
Water standard state calculator.
Definition PDSS_HKFT.h:213
double m_Mu0_tr_pr
Value of the Absolute Gibbs Free Energy NIST scale at T_r and P_r.
Definition PDSS_HKFT.h:248
void set_c(double *c)
Set "c" coefficients (array of 2 elements).
double m_a2
Input a2 coefficient (cal gmol-1)
Definition PDSS_HKFT.h:260
void setDeltaH0(double dh0)
Set enthalpy of formation at Pr, Tr [J/kmol].
double cp_R_ref() const override
Return the molar heat capacity divided by R at reference pressure.
double cp_R() const override
Return the molar const pressure heat capacity divided by RT.
Definition PDSS.cpp:179
double entropy_R() const override
Return the standard state entropy divided by RT.
Definition PDSS.cpp:169
double enthalpy_RT() const override
Return the standard state molar enthalpy divided by RT.
Definition PDSS.cpp:164
double gibbs_RT() const override
Return the molar Gibbs free energy divided by RT.
Definition PDSS.cpp:174
Class for the liquid water pressure dependent standard state.
Definition PDSS_Water.h:50
double thermalExpansionCoeff() const override
Return the volumetric thermal expansion coefficient. Units: 1/K.
double isothermalCompressibility() const
Returns the isothermal compressibility. Units: 1/Pa.
double dthermalExpansionCoeffdT() const
Return the derivative of the volumetric thermal expansion coefficient.
void setState_TP(double temp, double pres) override
Set the internal temperature and pressure.
double pref_safe(double temp) const
Returns a reference pressure value that can be safely calculated by the underlying real equation of s...
double density() const override
Return the standard state density at standard state.
virtual void setTemperature(double temp)
Set the internal temperature.
Definition PDSS.cpp:138
virtual void initThermo()
Initialization routine.
Definition PDSS.h:383
double m_temp
Current temperature used by the PDSS object.
Definition PDSS.h:398
double m_pres
State of the system - pressure.
Definition PDSS.h:401
double m_mw
Molecular Weight of the species.
Definition PDSS.h:413
AnyMap m_input
Input data supplied via setParameters.
Definition PDSS.h:417
virtual void getParameters(AnyMap &eosNode) const
Store the parameters needed to reconstruct a copy of this PDSS object.
Definition PDSS.h:392
virtual void setPressure(double pres)
Sets the pressure in the object.
Definition PDSS.cpp:128
size_t elementIndex(const string &name) const
Return the index of element named 'name'.
Definition Phase.cpp:55
string speciesName(size_t k) const
Name of the species with index k.
Definition Phase.cpp:142
double nAtoms(size_t k, size_t m) const
Number of atoms of element m in species k.
Definition Phase.cpp:103
size_t nElements() const
Number of elements.
Definition Phase.cpp:30
string elementName(size_t m) const
Name of the element with index m.
Definition Phase.cpp:49
double charge(size_t k) const
Dimensionless electrical charge of a single molecule of species k The charge is normalized by the the...
Definition Phase.h:538
double entropyElement298(size_t m) const
Entropy of the element in its standard state at 298 K and 1 bar.
Definition Phase.cpp:75
double convertFrom(double value, const string &src) const
Convert value from the specified src units to units appropriate for this unit system (defined by setD...
Definition Units.cpp:570
double convertTo(double value, const string &dest) const
Convert value to the specified dest units from the appropriate units for this unit system (defined by...
Definition Units.cpp:554
void setDefaults(std::initializer_list< string > units)
Set the default units to convert from when explicit units are not provided.
Definition Units.cpp:428
This file contains definitions for utility functions and text for modules, inputfiles and logging,...
const double OneAtm
One atmosphere [Pa].
Definition ct_defs.h:96
void warn_user(const string &method, const string &msg, const Args &... args)
Print a user warning raised from method as CanteraWarning.
Definition global.h:267
Namespace for the Cantera kernel.
Definition AnyMap.cpp:564
const size_t npos
index returned by functions to indicate "no position"
Definition ct_defs.h:180
Contains declarations for string manipulation functions within Cantera.