1 #include "cantera/thermo/ThermoFactory.h"
2 #include "cantera/thermo/WaterSSTP.h"
3 #include <iostream>
4
5 using namespace std;
6 using namespace Cantera;
7
tvalue(double val,double atol=1.0E-9)8 double tvalue(double val, double atol = 1.0E-9)
9 {
10 double rval = val;
11 if (fabs(val) < atol) {
12 rval = 0.0;
13 }
14 return rval;
15 }
16
main()17 int main()
18 {
19 #if defined(_MSC_VER) && _MSC_VER < 1900
20 _set_output_format(_TWO_DIGIT_EXPONENT);
21 #endif
22 double pres;
23 try {
24 ThermoPhase* w = newPhase("liquidvapor.yaml", "liquid-water-IAPWS95");
25 (dynamic_cast<WaterSSTP*>(w))->_allowGasPhase(true);
26
27 /*
28 * Print out the triple point conditions
29 */
30 double temp = 273.16;
31 pres = w->satPressure(temp);
32 printf("psat(%g) = %g\n", temp, pres);
33 double presLow = 1.0E-2;
34 temp = 298.15;
35 double oneBar = 1.0E5;
36 double vol;
37
38 printf("Comparisons to NIST: (see http://webbook.nist.gov):\n\n");
39
40 w->setDensity(1.0E-8);
41 w->setState_TP(temp, presLow);
42 double h = w->enthalpy_mole();
43 printf("H0(298.15) = %g J/kmol\n", h);
44 double h298 = h;
45
46 double s = w->entropy_mole();
47 s -= GasConstant * log(oneBar/presLow);
48 printf("S0(298.15) = %g J/kmolK\n", s);
49
50
51 double T[20];
52 T[0] = 298.15;
53 T[1] = 500.;
54 T[2] = 600.;
55 T[3] = 1000.;
56
57 double Cp0, delh0, delg0, g;
58 double Cp0_ss;
59 printf("\nIdeal Gas Standard State:\n");
60 printf(" T Cp0 S0 "
61 " -(G0-H298)/T H0-H298\n");
62 printf(" (K) (J/molK) (J/molK) "
63 " (J/molK) (kJ/mol)\n");
64 for (int i = 0; i < 4; i++) {
65 temp = T[i];
66 w->setState_TP(temp, presLow);
67 h = w->enthalpy_mole();
68 delh0 = tvalue(h - h298, 1.0E-6);
69 g = w->gibbs_mole();
70 delg0 = (g - h298)/temp + GasConstant * log(oneBar/presLow);
71 Cp0 = w->cp_mole();
72 {
73 w->getCp_R(&Cp0_ss);
74 Cp0_ss *= GasConstant;
75 if (fabs(Cp0_ss - Cp0) > 1.0E-5) {
76 printf("Inconsistency!\n");
77 exit(-1);
78 }
79 }
80 s = w->entropy_mole();
81 s -= GasConstant * log(oneBar/presLow);
82 printf("%10g %10g %13.4f %13.4f %13.4f\n", temp, Cp0*1.0E-3, s*1.0E-3,
83 -delg0*1.0E-3, delh0*1.0E-6);
84 }
85 printf("\n\n");
86
87 temp = 298.15;
88 w->setDensity(1000.);
89 w->setState_TP(temp, oneBar);
90 h = w->enthalpy_mole();
91 printf("H_liq(298.15, onebar) = %g J/kmol\n", h);
92 double h298l = h;
93 s = w->entropy_mole();
94 printf("S_liq(298.15, onebar) = %g J/kmolK\n", s);
95
96
97 T[0] = 273.19;
98 T[1] = 298.15;
99 T[2] = 300.;
100 T[3] = 373.15;
101 T[4] = 400.;
102 T[5] = 500.;
103 printf("\nLiquid 1bar or psat Standard State\n");
104 printf(" T press psat Cp0 S0 "
105 " -(G0-H298)/T H0-H298\n");
106 printf(" (K) (bar) (bar) (J/molK) (J/molK)"
107 " (J/molK) (kJ/mol)\n");
108
109 for (int i = 0; i < 6; i++) {
110 temp = T[i];
111 double psat = w->satPressure(temp);
112 double press = oneBar;
113 if (psat > press) {
114 press = psat*1.002;
115 }
116 w->setState_TP(temp, press);
117 h = w->enthalpy_mole();
118 delh0 = tvalue(h - h298l, 1.0E-6);
119 g = w->gibbs_mole();
120 delg0 = (g - h298l)/temp;
121 Cp0 = w->cp_mole();
122 s = w->entropy_mole();
123 printf("%10g %10g %12g %13.4f %13.4f %13.4f %13.4f\n", temp, press*1.0E-5,
124 psat*1.0E-5,
125 Cp0*1.0E-3, s*1.0E-3,
126 -delg0*1.0E-3, delh0*1.0E-6);
127 }
128
129 printf("\nLiquid Densities:\n");
130 printf(" T press psat Density molarVol "
131 "\n");
132 printf(" (K) (bar) (bar) (kg/m3) (m3/kmol)"
133 "\n");
134 for (int i = 0; i < 6; i++) {
135 temp = T[i];
136 double psat = w->satPressure(temp);
137 double press = oneBar;
138 if (psat > press) {
139 press = psat*1.002;
140 }
141 w->setState_TP(temp, press);
142 double d = w->density();
143 double mw = w->molecularWeight(0);
144 double vbar = mw/d;
145 printf("%10g %10g %12g %13.4f %13.4f\n", temp, press*1.0E-5,
146 psat*1.0E-5, d, vbar);
147
148 }
149
150
151 printf("\nLiquid 1bar or psat State: Partial Molar Quantities\n");
152 printf(" T press psat Cpbar Sbar "
153 " -(G0-H298)/T H0-H298 Volume\n");
154 printf(" (K) (bar) (bar) (J/molK) (J/molK)"
155 " (J/molK) (kJ/mol) m3/kmol\n");
156
157 for (int i = 0; i < 6; i++) {
158 temp = T[i];
159 double psat = w->satPressure(temp);
160 double press = oneBar;
161 if (psat > press) {
162 press = psat*1.002;
163 }
164 w->setState_TP(temp, press);
165 w->getPartialMolarEnthalpies(&h);
166 delh0 = tvalue(h - h298l, 1.0E-6);
167 w->getChemPotentials(&g);
168 delg0 = (g - h298l)/temp;
169 w->getPartialMolarCp(&Cp0);
170 w->getPartialMolarEntropies(&s);
171 w->getPartialMolarVolumes(&vol);
172 printf("%10g %10g %12g %13.4f %13.4f %13.4f %13.4f %13.4f\n", temp, press*1.0E-5,
173 psat*1.0E-5,
174 Cp0*1.0E-3, s*1.0E-3,
175 -delg0*1.0E-3, delh0*1.0E-6, vol);
176 }
177
178 printf("\nLiquid 1bar or psat State: Standard State Quantities\n");
179 printf(" T press psat Cpbar Sbar "
180 " -(G0-H298)/T H0-H298 Volume\n");
181 printf(" (K) (bar) (bar) (J/molK) (J/molK)"
182 " (J/molK) (kJ/mol) m3/kmol\n");
183
184 for (int i = 0; i < 6; i++) {
185 temp = T[i];
186 double psat = w->satPressure(temp);
187 double press = oneBar;
188 if (psat > press) {
189 press = psat*1.002;
190 }
191 w->setState_TP(temp, press);
192 w->getEnthalpy_RT(&h);
193 h *= temp * GasConstant;
194 delh0 = tvalue(h - h298l, 1.0E-6);
195 w->getStandardChemPotentials(&g);
196 delg0 = (g - h298l)/temp;
197 w->getCp_R(&Cp0);
198 Cp0 *= GasConstant;
199 w->getEntropy_R(&s);
200 s *= GasConstant;
201 w->getStandardVolumes(&vol);
202 printf("%10g %10g %12g %13.4f %13.4f %13.4f %13.4f %13.4f\n", temp, press*1.0E-5,
203 psat*1.0E-5,
204 Cp0*1.0E-3, s*1.0E-3,
205 -delg0*1.0E-3, delh0*1.0E-6, vol);
206 }
207
208 printf("\nLiquid 1bar or psat State: Reference State Quantities (Always 1 atm no matter what system pressure is)\n");
209 printf(" T press psat Cpbar Sbar "
210 " -(G0-H298)/T H0-H298 Volume\n");
211 printf(" (K) (bar) (bar) (J/molK) (J/molK)"
212 " (J/molK) (kJ/mol) m3/kmol\n");
213
214 for (int i = 0; i < 6; i++) {
215 temp = T[i];
216 double psat = w->satPressure(temp);
217 double press = oneBar;
218 if (psat > press) {
219 press = psat*1.002;
220 }
221 w->setState_TP(temp, press);
222 w->getEnthalpy_RT_ref(&h);
223 h *= temp * GasConstant;
224 delh0 = tvalue(h - h298l, 1.0E-6);
225 w->getGibbs_ref(&g);
226 delg0 = (g - h298l)/temp;
227 w->getCp_R_ref(&Cp0);
228 Cp0 *= GasConstant;
229 w->getEntropy_R_ref(&s);
230 s *= GasConstant;
231 w->getStandardVolumes_ref(&vol);
232 printf("%10g %10g %12g %13.4f %13.4f %13.4f %13.4f %13.4f\n", temp, press*1.0E-5,
233 psat*1.0E-5,
234 Cp0*1.0E-3, s*1.0E-3,
235 -delg0*1.0E-3, delh0*1.0E-6, vol);
236 }
237
238 printf("\nLiquid 1 atm: Standard State Quantities - Should agree with table above\n");
239 printf(" T press psat Cpbar Sbar "
240 " -(G0-H298)/T H0-H298 Volume\n");
241 printf(" (K) (bar) (bar) (J/molK) (J/molK)"
242 " (J/molK) (kJ/mol) m3/kmol\n");
243
244 for (int i = 0; i < 6; i++) {
245 temp = T[i];
246 double psat = w->satPressure(temp);
247 double press = OneAtm;
248 w->setState_TP(temp, press);
249 w->getEnthalpy_RT(&h);
250 h *= temp * GasConstant;
251 delh0 = tvalue(h - h298l, 1.0E-6);
252 w->getStandardChemPotentials(&g);
253 delg0 = (g - h298l)/temp;
254 w->getCp_R(&Cp0);
255 Cp0 *= GasConstant;
256 w->getEntropy_R(&s);
257 s *= GasConstant;
258 w->getStandardVolumes(&vol);
259 printf("%10g %10g %12g %13.4f %13.4f %13.4f %13.4f %13.4f\n", temp, press*1.0E-5,
260 psat*1.0E-5,
261 Cp0*1.0E-3, s*1.0E-3,
262 -delg0*1.0E-3, delh0*1.0E-6, vol);
263 }
264 printf("\n\nTable of increasing Enthalpy at 1 atm\n\n");
265 double dens;
266 printf(" Enthalpy, Temperature, x_Vapor, Density, Entropy_mass, Gibbs_mass\n");
267 w->setState_TP(298., OneAtm);
268 double Hset = w->enthalpy_mass();
269 double vapFrac = w->vaporFraction();
270 double Tcalc = w->temperature();
271 double Scalc = w->entropy_mass();
272 double Gcalc = w->gibbs_mass();
273 dens = w->density();
274 printf(" %10g, %10g, %10g, %11.5g, %11.5g, %11.5g\n", Hset , Tcalc, vapFrac, dens, Scalc, Gcalc);
275 w->setState_HP(Hset, OneAtm);
276 vapFrac = w->vaporFraction();
277 Tcalc = w->temperature();
278 dens = w->density();
279 Scalc = w->entropy_mass();
280 Gcalc = w->gibbs_mass();
281 printf(" %10g, %10g, %10g, %11.5g, %11.5g, %11.5g\n", Hset , Tcalc, vapFrac, dens, Scalc, Gcalc);
282
283 double deltaH = 100000.;
284 for (int i = 0; i < 40; i++) {
285 Hset += deltaH;
286 try {
287 w->setState_HP(Hset, OneAtm);
288 } catch (CanteraError&) {
289 printf(" %10g, -> Failed to converge, beyond the spinodal probably \n\n", Hset);
290 break;
291 }
292 vapFrac = w->vaporFraction();
293 Tcalc = w->temperature();
294 dens = w->density();
295 Scalc = w->entropy_mass();
296 Gcalc = w->gibbs_mass();
297 printf(" %10g, %10g, %10g, %11.5g, %11.5g, %11.5g\n", Hset , Tcalc, vapFrac, dens, Scalc, Gcalc);
298 }
299
300
301
302 printf("Critical Temp = %10.3g K\n", w->critTemperature());
303 printf("Critical Pressure = %10.3g atm\n", w->critPressure()/OneAtm);
304 printf("Critical Dens = %10.3g kg/m3\n", w->critDensity());
305
306 delete w;
307 } catch (CanteraError& err) {
308 std::cout << err.what() << std::endl;
309 Cantera::appdelete();
310 return -1;
311 }
312
313
314 return 0;
315 }
316