PDSS_HKFT.cpp Source File#

Cantera: PDSS_HKFT.cpp Source File
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
12#include "cantera/thermo/PDSS_HKFT.h"
13#include "cantera/thermo/PDSS_Water.h"
14#include "cantera/thermo/VPStandardStateTP.h"
15#include "cantera/base/stringUtils.h"
16#include "cantera/base/global.h"
17
18namespace Cantera
19{
20// Set the default to error exit if there is an input file inconsistency
21int PDSS_HKFT::s_InputInconsistencyErrorExit = 1;
22
23PDSS_HKFT::PDSS_HKFT()
24{
25 m_pres = OneAtm;
26 m_units.setDefaults({"cm", "gmol", "cal", "bar"});
27}
28
29double PDSS_HKFT::enthalpy_mole() const
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
36double PDSS_HKFT::intEnergy_mole() const
37{
38 return enthalpy_RT() - molarVolume() * m_pres;
39}
40
41double PDSS_HKFT::entropy_mole() const
42{
43 return m_units.convertTo(m_Entrop_tr_pr, "J/kmol") + deltaS();
44}
45
46double PDSS_HKFT::gibbs_mole() const
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
107double PDSS_HKFT::molarVolume() const
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
155double PDSS_HKFT::gibbs_RT_ref() const
156{
157 double m_psave = m_pres;
158 m_pres = m_waterSS->pref_safe(m_temp);
159 double ee = gibbs_RT();
160 m_pres = m_psave;
161 return ee;
162}
163
164double PDSS_HKFT::enthalpy_RT_ref() const
165{
166 double m_psave = m_pres;
167 m_pres = m_waterSS->pref_safe(m_temp);
168 double hh = enthalpy_RT();
169 m_pres = m_psave;
170 return hh;
171}
172
173double PDSS_HKFT::entropy_R_ref() const
174{
175 double m_psave = m_pres;
176 m_pres = m_waterSS->pref_safe(m_temp);
177 double ee = entropy_R();
178 m_pres = m_psave;
179 return ee;
180}
181
182double PDSS_HKFT::cp_R_ref() const
183{
184 double m_psave = m_pres;
185 m_pres = m_waterSS->pref_safe(m_temp);
186 double ee = cp_R();
187 m_pres = m_psave;
188 return ee;
189}
190
191double PDSS_HKFT::molarVolume_ref() const
192{
193 double m_psave = m_pres;
194 m_pres = m_waterSS->pref_safe(m_temp);
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
206void PDSS_HKFT::initThermo()
207{
208 PDSS::initThermo();
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.
237 m_charge_j = m_tp->charge(m_spindex);
238 if (std::isnan(m_deltaH_formation_tr_pr)) {
239 convertDGFormation();
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");
245 m_deltaG_formation_tr_pr = m_units.convertFrom(m_Mu0_tr_pr, "J/kmol");
246 double tmp = m_Mu0_tr_pr;
247 convertDGFormation();
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)) {
252 convertDGFormation();
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);
263 m_waterSS->setState_TP(m_temp, m_pres);
264 m_densWaterSS = m_waterSS->density();
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;
271 convertDGFormation();
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);
282 if (s_InputInconsistencyErrorExit) {
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
343void PDSS_HKFT::getParameters(AnyMap& eosNode) const
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);
505 m_densWaterSS = m_waterSS->density();
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
575void PDSS_HKFT::convertDGFormation()
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}
ENTROPY298_UNKNOWN
#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
PDSS_HKFT.h
Declarations for the class PDSS_HKFT (pressure dependent standard state) which handles calculations f...
PDSS_Water.h
Implementation of a pressure dependent standard state virtual function for a Pure Water Phase (see Sp...
VPStandardStateTP.h
Header file for a derived class of ThermoPhase that handles variable pressure standard state methods ...
Cantera::AnyMap
A map of string keys to values whose type can vary at runtime.
Definition AnyMap.h:427
Cantera::AnyMap::hasKey
bool hasKey(const string &key) const
Returns true if the map contains an item named key.
Definition AnyMap.cpp:1423
Cantera::AnyMap::units
const UnitSystem & units() const
Return the default units that should be used to convert stored values.
Definition AnyMap.h:630
Cantera::AnyMap::convert
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
Cantera::AnyValue
A wrapper for a variable whose type is determined at runtime.
Definition AnyMap.h:86
Cantera::CanteraError
Base class for exceptions thrown by Cantera classes.
Definition ctexceptions.h:66
Cantera::NotImplementedError
An error indicating that an unimplemented function has been called.
Definition ctexceptions.h:195
Cantera::PDSS_HKFT::setOmega
void setOmega(double omega)
Set omega [J/kmol].
Definition PDSS_HKFT.cpp:339
Cantera::PDSS_HKFT::molarVolume
double molarVolume() const override
Return the molar volume at standard state.
Definition PDSS_HKFT.cpp:107
Cantera::PDSS_HKFT::s_InputInconsistencyErrorExit
static int s_InputInconsistencyErrorExit
Static variable determining error exiting.
Definition PDSS_HKFT.h:297
Cantera::PDSS_HKFT::m_Y_pr_tr
double m_Y_pr_tr
y = dZdT = 1/(esp*esp) desp/dT at 298.15 and 1 bar
Definition PDSS_HKFT.h:278
Cantera::PDSS_HKFT::m_a4
double m_a4
Input a4 coefficient (cal K gmol-1)
Definition PDSS_HKFT.h:266
Cantera::PDSS_HKFT::enthalpy_mole
double enthalpy_mole() const override
Return the molar enthalpy in units of J kmol-1.
Definition PDSS_HKFT.cpp:29
Cantera::PDSS_HKFT::m_waterProps
unique_ptr< WaterProps > m_waterProps
Pointer to the water property calculator.
Definition PDSS_HKFT.h:221
Cantera::PDSS_HKFT::deltaS
double deltaS() const
Main routine that actually calculates the entropy difference between the reference state at Tr,...
Definition PDSS_HKFT.cpp:398
Cantera::PDSS_HKFT::f
double f(const double temp, const double pres, const int ifunc=0) const
Difference function f appearing in the formulation.
Definition PDSS_HKFT.cpp:466
Cantera::PDSS_HKFT::m_spindex
size_t m_spindex
Index of this species within the parent phase.
Definition PDSS_HKFT.h:101
Cantera::PDSS_HKFT::m_densWaterSS
double m_densWaterSS
density of standard-state water. internal temporary variable
Definition PDSS_HKFT.h:218
Cantera::PDSS_HKFT::LookupGe
double LookupGe(const string &elemName)
Function to look up Element Free Energies.
Definition PDSS_HKFT.cpp:561
Cantera::PDSS_HKFT::gibbs_RT_ref
double gibbs_RT_ref() const override
Return the molar Gibbs free energy divided by RT at reference pressure.
Definition PDSS_HKFT.cpp:155
Cantera::PDSS_HKFT::m_tp
VPStandardStateTP * m_tp
Parent VPStandardStateTP (ThermoPhase) object.
Definition PDSS_HKFT.h:100
Cantera::PDSS_HKFT::m_Entrop_tr_pr
double m_Entrop_tr_pr
Input value of S_j at Tr and Pr (cal gmol-1 K-1)
Definition PDSS_HKFT.h:254
Cantera::PDSS_HKFT::m_deltaG_formation_tr_pr
double m_deltaG_formation_tr_pr
Input value of deltaG of Formation at Tr and Pr (cal gmol-1)
Definition PDSS_HKFT.h:230
Cantera::PDSS_HKFT::setDeltaG0
void setDeltaG0(double dg0)
Set Gibbs free energy of formation at Pr, Tr [J/kmol].
Definition PDSS_HKFT.cpp:319
Cantera::PDSS_HKFT::initThermo
void initThermo() override
Initialization routine.
Definition PDSS_HKFT.cpp:206
Cantera::PDSS_HKFT::m_a3
double m_a3
Input a3 coefficient (cal K gmol-1 bar-1)
Definition PDSS_HKFT.h:263
Cantera::PDSS_HKFT::gstar
double gstar(const double temp, const double pres, const int ifunc=0) const
Evaluate the Gstar value appearing in the HKFT formulation.
Definition PDSS_HKFT.cpp:553
Cantera::PDSS_HKFT::m_charge_j
double m_charge_j
Charge of the ion.
Definition PDSS_HKFT.h:290
Cantera::PDSS_HKFT::PDSS_HKFT
PDSS_HKFT()
Default Constructor.
Definition PDSS_HKFT.cpp:23
Cantera::PDSS_HKFT::entropy_R_ref
double entropy_R_ref() const override
Return the molar entropy divided by R at reference pressure.
Definition PDSS_HKFT.cpp:173
Cantera::PDSS_HKFT::deltaG
double deltaG() const
Main routine that actually calculates the Gibbs free energy difference between the reference state at...
Definition PDSS_HKFT.cpp:365
Cantera::PDSS_HKFT::enthalpy_RT_ref
double enthalpy_RT_ref() const override
Return the molar enthalpy divided by RT at reference pressure.
Definition PDSS_HKFT.cpp:164
Cantera::PDSS_HKFT::setS0
void setS0(double s0)
Set entropy of formation at Pr, Tr [J/kmol/K].
Definition PDSS_HKFT.cpp:323
Cantera::PDSS_HKFT::bg
double bg(const double temp, const int ifunc=0) const
Internal formula for the calculation of b_g()
Definition PDSS_HKFT.cpp:452
Cantera::PDSS_HKFT::getParameters
void getParameters(AnyMap &eosNode) const override
Store the parameters needed to reconstruct a copy of this PDSS object.
Definition PDSS_HKFT.cpp:343
Cantera::PDSS_HKFT::m_c2
double m_c2
Input c2 coefficient (cal K gmol-1)
Definition PDSS_HKFT.h:272
Cantera::PDSS_HKFT::m_presR_bar
double m_presR_bar
Reference pressure is 1 atm in units of bar= 1.0132.
Definition PDSS_HKFT.h:284
Cantera::PDSS_HKFT::set_a
void set_a(double *a)
Set "a" coefficients (array of 4 elements).
Definition PDSS_HKFT.cpp:327
Cantera::PDSS_HKFT::intEnergy_mole
double intEnergy_mole() const override
Return the molar internal Energy in units of J kmol-1.
Definition PDSS_HKFT.cpp:36
Cantera::PDSS_HKFT::m_domega_jdT_prtr
double m_domega_jdT_prtr
small value that is not quite zero
Definition PDSS_HKFT.h:287
Cantera::PDSS_HKFT::entropy_mole
double entropy_mole() const override
Return the molar entropy in units of J kmol-1 K-1.
Definition PDSS_HKFT.cpp:41
Cantera::PDSS_HKFT::m_omega_pr_tr
double m_omega_pr_tr
Input omega_pr_tr coefficient(cal gmol-1)
Definition PDSS_HKFT.h:275
Cantera::PDSS_HKFT::g
double g(const double temp, const double pres, const int ifunc=0) const
function g appearing in the formulation
Definition PDSS_HKFT.cpp:500
Cantera::PDSS_HKFT::setState_TP
void setState_TP(double temp, double pres) override
Set the internal temperature and pressure.
Definition PDSS_HKFT.cpp:200
Cantera::PDSS_HKFT::m_Z_pr_tr
double m_Z_pr_tr
Z = -1 / relEpsilon at 298.15 and 1 bar.
Definition PDSS_HKFT.h:281
Cantera::PDSS_HKFT::convertDGFormation
void convertDGFormation()
Translate a Gibbs free energy of formation value to a NIST-based Chemical potential.
Definition PDSS_HKFT.cpp:575
Cantera::PDSS_HKFT::cp_mole
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
Cantera::PDSS_HKFT::m_deltaH_formation_tr_pr
double m_deltaH_formation_tr_pr
Input value of deltaH of Formation at Tr and Pr (cal gmol-1)
Definition PDSS_HKFT.h:239
Cantera::PDSS_HKFT::m_a1
double m_a1
Input a1 coefficient (cal gmol-1 bar-1)
Definition PDSS_HKFT.h:257
Cantera::PDSS_HKFT::ag
double ag(const double temp, const int ifunc=0) const
Internal formula for the calculation of a_g()
Definition PDSS_HKFT.cpp:438
Cantera::PDSS_HKFT::density
double density() const override
Return the standard state density at standard state.
Definition PDSS_HKFT.cpp:150
Cantera::PDSS_HKFT::gibbs_mole
double gibbs_mole() const override
Return the molar Gibbs free energy in units of J kmol-1.
Definition PDSS_HKFT.cpp:46
Cantera::PDSS_HKFT::m_c1
double m_c1
Input c1 coefficient (cal gmol-1 K-1)
Definition PDSS_HKFT.h:269
Cantera::PDSS_HKFT::molarVolume_ref
double molarVolume_ref() const override
Return the molar volume at reference pressure.
Definition PDSS_HKFT.cpp:191
Cantera::PDSS_HKFT::m_waterSS
PDSS_Water * m_waterSS
Water standard state calculator.
Definition PDSS_HKFT.h:213
Cantera::PDSS_HKFT::m_Mu0_tr_pr
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
Cantera::PDSS_HKFT::set_c
void set_c(double *c)
Set "c" coefficients (array of 2 elements).
Definition PDSS_HKFT.cpp:334
Cantera::PDSS_HKFT::m_a2
double m_a2
Input a2 coefficient (cal gmol-1)
Definition PDSS_HKFT.h:260
Cantera::PDSS_HKFT::setDeltaH0
void setDeltaH0(double dh0)
Set enthalpy of formation at Pr, Tr [J/kmol].
Definition PDSS_HKFT.cpp:315
Cantera::PDSS_HKFT::cp_R_ref
double cp_R_ref() const override
Return the molar heat capacity divided by R at reference pressure.
Definition PDSS_HKFT.cpp:182
Cantera::PDSS_Molar::cp_R
double cp_R() const override
Return the molar const pressure heat capacity divided by RT.
Definition PDSS.cpp:179
Cantera::PDSS_Molar::entropy_R
double entropy_R() const override
Return the standard state entropy divided by RT.
Definition PDSS.cpp:169
Cantera::PDSS_Molar::enthalpy_RT
double enthalpy_RT() const override
Return the standard state molar enthalpy divided by RT.
Definition PDSS.cpp:164
Cantera::PDSS_Molar::gibbs_RT
double gibbs_RT() const override
Return the molar Gibbs free energy divided by RT.
Definition PDSS.cpp:174
Cantera::PDSS_Water
Class for the liquid water pressure dependent standard state.
Definition PDSS_Water.h:50
Cantera::PDSS_Water::thermalExpansionCoeff
double thermalExpansionCoeff() const override
Return the volumetric thermal expansion coefficient. Units: 1/K.
Definition PDSS_Water.cpp:160
Cantera::PDSS_Water::isothermalCompressibility
double isothermalCompressibility() const
Returns the isothermal compressibility. Units: 1/Pa.
Definition PDSS_Water.cpp:181
Cantera::PDSS_Water::dthermalExpansionCoeffdT
double dthermalExpansionCoeffdT() const
Return the derivative of the volumetric thermal expansion coefficient.
Definition PDSS_Water.cpp:165
Cantera::PDSS_Water::setState_TP
void setState_TP(double temp, double pres) override
Set the internal temperature and pressure.
Definition PDSS_Water.cpp:218
Cantera::PDSS_Water::pref_safe
double pref_safe(double temp) const
Returns a reference pressure value that can be safely calculated by the underlying real equation of s...
Definition PDSS_Water.cpp:224
Cantera::PDSS_Water::density
double density() const override
Return the standard state density at standard state.
Definition PDSS_Water.cpp:207
Cantera::PDSS::setTemperature
virtual void setTemperature(double temp)
Set the internal temperature.
Definition PDSS.cpp:138
Cantera::PDSS::initThermo
virtual void initThermo()
Initialization routine.
Definition PDSS.h:383
Cantera::PDSS::m_temp
double m_temp
Current temperature used by the PDSS object.
Definition PDSS.h:398
Cantera::PDSS::m_pres
double m_pres
State of the system - pressure.
Definition PDSS.h:401
Cantera::PDSS::m_mw
double m_mw
Molecular Weight of the species.
Definition PDSS.h:413
Cantera::PDSS::m_input
AnyMap m_input
Input data supplied via setParameters.
Definition PDSS.h:417
Cantera::PDSS::getParameters
virtual void getParameters(AnyMap &eosNode) const
Store the parameters needed to reconstruct a copy of this PDSS object.
Definition PDSS.h:392
Cantera::PDSS::setPressure
virtual void setPressure(double pres)
Sets the pressure in the object.
Definition PDSS.cpp:128
Cantera::Phase::elementIndex
size_t elementIndex(const string &name) const
Return the index of element named 'name'.
Definition Phase.cpp:55
Cantera::Phase::speciesName
string speciesName(size_t k) const
Name of the species with index k.
Definition Phase.cpp:142
Cantera::Phase::nAtoms
double nAtoms(size_t k, size_t m) const
Number of atoms of element m in species k.
Definition Phase.cpp:103
Cantera::Phase::nElements
size_t nElements() const
Number of elements.
Definition Phase.cpp:30
Cantera::Phase::elementName
string elementName(size_t m) const
Name of the element with index m.
Definition Phase.cpp:49
Cantera::Phase::charge
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
Cantera::Phase::entropyElement298
double entropyElement298(size_t m) const
Entropy of the element in its standard state at 298 K and 1 bar.
Definition Phase.cpp:75
Cantera::UnitSystem::convertFrom
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
Cantera::UnitSystem::convertTo
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
Cantera::UnitSystem::setDefaults
void setDefaults(std::initializer_list< string > units)
Set the default units to convert from when explicit units are not provided.
Definition Units.cpp:428
global.h
This file contains definitions for utility functions and text for modules, inputfiles and logging,...
Cantera::OneAtm
const double OneAtm
One atmosphere [Pa].
Definition ct_defs.h:96
Cantera::warn_user
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
Cantera
Namespace for the Cantera kernel.
Definition AnyMap.cpp:564
Cantera::npos
const size_t npos
index returned by functions to indicate "no position"
Definition ct_defs.h:180
stringUtils.h
Contains declarations for string manipulation functions within Cantera.
Generated by doxygen 1.9.7