Cantera  3.1.0a1
Loading...
Searching...
No Matches
WaterPropsIAPWS.cpp
Go to the documentation of this file.
1/**
2 * @file WaterPropsIAPWS.cpp
3 * Definitions for a class for calculating the equation of state of water
4 * from the IAPWS 1995 Formulation based on the steam tables thermodynamic
5 * basis (See class @link Cantera::WaterPropsIAPWS WaterPropsIAPWS@endlink).
6 */
7
8// This file is part of Cantera. See License.txt in the top-level directory or
9// at https://cantera.org/license.txt for license and copyright information.
10
14#include "cantera/base/global.h"
15
16namespace Cantera
17{
18// Critical Point values of water in mks units
19
20//! Critical Temperature value (kelvin)
21const double T_c = 647.096;
22//! Critical Pressure (Pascals)
23static const double P_c = 22.064E6;
24//! Value of the Density at the critical point (kg m-3)
25const double Rho_c = 322.;
26
27static const double R_water = 461.51805; // J/kg/K (Eq. 6.3)
28
29void WaterPropsIAPWS::calcDim(double temperature, double rho)
30{
32 delta = rho / Rho_c;
33
34 // Determine the internal state
35 if (temperature > T_c) {
36 iState = WATER_SUPERCRIT;
37 } else {
38 if (delta < 1.0) {
39 iState = WATER_GAS;
40 } else {
41 iState = WATER_LIQUID;
42 }
43 }
44}
45
47{
48 double retn = m_phi.pressureM_rhoRT(tau, delta);
49 double rho = delta * Rho_c;
50 double temperature = T_c / tau;
51 return retn * rho * R_water * temperature;
52}
53
54double WaterPropsIAPWS::density(double temperature, double pressure,
55 int phase, double rhoguess)
56{
57 if (fabs(pressure - P_c) / P_c < 1.e-8 &&
58 fabs(temperature - T_c) / T_c < 1.e-8) {
59 // Catch critical point, as no solution is found otherwise
61 return Rho_c;
62 }
63 double deltaGuess = 0.0;
64 if (rhoguess == -1.0) {
65 if (phase != -1) {
66 if (temperature > T_c) {
67 rhoguess = pressure / (R_water * temperature);
68 } else {
69 if (phase == WATER_GAS || phase == WATER_SUPERCRIT) {
70 rhoguess = pressure / (R_water * temperature);
71 } else if (phase == WATER_LIQUID) {
72 // Provide a guess about the liquid density that is
73 // relatively high -> convergence from above seems robust.
74 rhoguess = 1000.;
75 } else if (phase == WATER_UNSTABLELIQUID || phase == WATER_UNSTABLEGAS) {
76 throw CanteraError("WaterPropsIAPWS::density",
77 "Unstable Branch finder is untested");
78 } else {
79 throw CanteraError("WaterPropsIAPWS::density",
80 "unknown state: {}", phase);
81 }
82 }
83 } else {
84 // Assume the Gas phase initial guess, if nothing is specified to
85 // the routine
86 rhoguess = pressure / (R_water * temperature);
87 }
88 }
89 double p_red = pressure / (R_water * temperature * Rho_c);
90 deltaGuess = rhoguess / Rho_c;
91 setState_TD(temperature, rhoguess);
92 double delta_retn = m_phi.dfind(p_red, tau, deltaGuess);
93 if (delta_retn <= 0) {
94 // No solution found for first initial guess; perturb initial guess once
95 // to avoid spurious failures (band-aid fix)
96 delta_retn = m_phi.dfind(p_red, tau, 0.9 * deltaGuess);
97 }
98 double density_retn;
99 if (delta_retn > 0.0) {
100 delta = delta_retn;
101
102 // Dimensionalize the density before returning
103 density_retn = delta_retn * Rho_c;
104
105 // Set the internal state -> this may be a duplication. However, let's
106 // just be sure.
107 setState_TD(temperature, density_retn);
108 } else {
109 density_retn = -1.0;
110 }
111 return density_retn;
112}
113
114double WaterPropsIAPWS::density_const(double pressure, int phase, double rhoguess) const
115{
116 double temperature = T_c / tau;
117 double deltaGuess = 0.0;
118 double deltaSave = delta;
119 if (rhoguess == -1.0) {
120 if (phase != -1) {
121 if (temperature > T_c) {
122 rhoguess = pressure / (R_water * temperature);
123 } else {
124 if (phase == WATER_GAS || phase == WATER_SUPERCRIT) {
125 rhoguess = pressure / (R_water * temperature);
126 } else if (phase == WATER_LIQUID) {
127 // Provide a guess about the liquid density that is
128 // relatively high -> convergence from above seems robust.
129 rhoguess = 1000.;
130 } else if (phase == WATER_UNSTABLELIQUID || phase == WATER_UNSTABLEGAS) {
131 throw CanteraError("WaterPropsIAPWS::density_const",
132 "Unstable Branch finder is untested");
133 } else {
134 throw CanteraError("WaterPropsIAPWS::density_const",
135 "unknown state: {}", phase);
136 }
137 }
138 } else {
139 // Assume the Gas phase initial guess, if nothing is specified to
140 // the routine
141 rhoguess = pressure / (R_water * temperature);
142 }
143 }
144 double p_red = pressure / (R_water * temperature * Rho_c);
145 deltaGuess = rhoguess / Rho_c;
146
147 delta = deltaGuess;
149
150 double delta_retn = m_phi.dfind(p_red, tau, deltaGuess);
151 double density_retn;
152 if (delta_retn > 0.0) {
153 delta = delta_retn;
154
155 // Dimensionalize the density before returning
156 density_retn = delta_retn * Rho_c;
157
158 } else {
159 density_retn = -1.0;
160 }
161
162 delta = deltaSave;
164 return density_retn;
165}
166
168{
169 return delta * Rho_c;
170}
171
173{
174 return T_c / tau;
175}
176
177double WaterPropsIAPWS::psat_est(double temperature) const
178{
179 // Formula and constants from: "NBS/NRC Steam Tables: Thermodynamic and
180 // Transport Properties and Computer Programs for Vapor and Liquid States of
181 // Water in SI Units". L. Haar, J. S. Gallagher, G. S. Kell. Hemisphere
182 // Publishing. 1984.
183 static const double A[8] = {
184 -7.8889166E0,
185 2.5514255E0,
186 -6.716169E0,
187 33.2239495E0,
188 -105.38479E0,
189 174.35319E0,
190 -148.39348E0,
191 48.631602E0
192 };
193 double ps;
194 if (temperature < 314.) {
195 double pl = 6.3573118E0 - 8858.843E0 / temperature
196 + 607.56335E0 * pow(temperature, -0.6);
197 ps = 0.1 * exp(pl);
198 } else {
199 double v = temperature / 647.25;
200 double w = fabs(1.0-v);
201 double b = 0.0;
202 for (int i = 0; i < 8; i++) {
203 double z = i + 1;
204 b += A[i] * pow(w, ((z+1.0)/2.0));
205 }
206 double q = b / v;
207 ps = 22.093*exp(q);
208 }
209
210 // Original correlation was in cgs. Convert to mks
211 ps *= 1.0E6;
212 return ps;
213}
214
216{
217 double dpdrho_val = dpdrho();
218 double dens = delta * Rho_c;
219 return 1.0 / (dens * dpdrho_val);
220}
221
223{
224 double retn = m_phi.dimdpdrho(tau, delta);
225 double temperature = T_c/tau;
226 return retn * R_water * temperature;
227}
228
230{
231 return m_phi.dimdpdT(tau, delta);
232}
233
235{
236 double kappa = isothermalCompressibility();
237 double beta = coeffPresExp();
238 double dens = delta * Rho_c;
239 return kappa * dens * R_water * beta;
240}
241
242void WaterPropsIAPWS::corr(double temperature, double pressure,
243 double& densLiq, double& densGas, double& delGRT)
244{
245 densLiq = density(temperature, pressure, WATER_LIQUID, densLiq);
246 if (densLiq <= 0.0) {
247 throw CanteraError("WaterPropsIAPWS::corr",
248 "Error occurred trying to find liquid density at (T,P) = {} {}",
250 }
251 setState_TD(temperature, densLiq);
252 double gibbsLiqRT = m_phi.gibbs_RT();
253
254 densGas = density(temperature, pressure, WATER_GAS, densGas);
255 if (densGas <= 0.0) {
256 throw CanteraError("WaterPropsIAPWS::corr",
257 "Error occurred trying to find gas density at (T,P) = {} {}",
259 }
260 setState_TD(temperature, densGas);
261 double gibbsGasRT = m_phi.gibbs_RT();
262
263 delGRT = gibbsLiqRT - gibbsGasRT;
264}
265
266void WaterPropsIAPWS::corr1(double temperature, double pressure,
267 double& densLiq, double& densGas, double& pcorr)
268{
269 densLiq = density(temperature, pressure, WATER_LIQUID, densLiq);
270 if (densLiq <= 0.0) {
271 throw CanteraError("WaterPropsIAPWS::corr1",
272 "Error occurred trying to find liquid density at (T,P) = {} {}",
274 }
275 setState_TD(temperature, densLiq);
276 double prL = m_phi.phiR();
277
278 densGas = density(temperature, pressure, WATER_GAS, densGas);
279 if (densGas <= 0.0) {
280 throw CanteraError("WaterPropsIAPWS::corr1",
281 "Error occurred trying to find gas density at (T,P) = {} {}",
283 }
284 setState_TD(temperature, densGas);
285 double prG = m_phi.phiR();
286 double rhs = (prL - prG) + log(densLiq/densGas);
287 rhs /= (1.0/densGas - 1.0/densLiq);
288 pcorr = rhs * R_water * temperature;
289}
290
291double WaterPropsIAPWS::psat(double temperature, int waterState)
292{
293 static int method = 1;
294 double densLiq = -1.0, densGas = -1.0, delGRT = 0.0;
295 double dp, pcorr;
296 if (temperature >= T_c) {
297 densGas = density(temperature, P_c, WATER_SUPERCRIT);
298 setState_TD(temperature, densGas);
299 return P_c;
300 }
301 double p = psat_est(temperature);
302 for (int i = 0; i < 30; i++) {
303 if (method == 1) {
304 corr(temperature, p, densLiq, densGas, delGRT);
305 double delV = 1.0/densLiq - 1.0/densGas;
306 dp = - delGRT * R_water * temperature / delV;
307 } else {
308 corr1(temperature, p, densLiq, densGas, pcorr);
309 dp = pcorr - p;
310 }
311 p += dp;
312
313 if ((method == 1) && delGRT < 1.0E-8) {
314 break;
315 } else {
316 if (fabs(dp/p) < 1.0E-9) {
317 break;
318 }
319 }
320 }
321 // Put the fluid in the desired end condition
322 if (waterState == WATER_LIQUID) {
323 setState_TD(temperature, densLiq);
324 } else if (waterState == WATER_GAS) {
325 setState_TD(temperature, densGas);
326 } else {
327 throw CanteraError("WaterPropsIAPWS::psat",
328 "unknown water state input: {}", waterState);
329 }
330 return p;
331}
332
333int WaterPropsIAPWS::phaseState(bool checkState) const
334{
335 if (checkState) {
336 if (tau <= 1.0) {
337 iState = WATER_SUPERCRIT;
338 } else {
339 double T = T_c / tau;
340 double rho = delta * Rho_c;
341 double rhoMidAtm = 0.5 * (OneAtm / (R_water * 373.15) + 1.0E3);
342 double rhoMid = Rho_c + (T - T_c) * (Rho_c - rhoMidAtm) / (T_c - 373.15);
343 int iStateGuess = WATER_LIQUID;
344 if (rho < rhoMid) {
345 iStateGuess = WATER_GAS;
346 }
347 double kappa = isothermalCompressibility();
348 if (kappa >= 0.0) {
349 iState = iStateGuess;
350 } else {
351 // When we are here we are between the spinodal curves
352 double rhoDel = rho * 1.000001;
353 double deltaSave = delta;
354 double deltaDel = rhoDel / Rho_c;
355 delta = deltaDel;
356 m_phi.tdpolycalc(tau, deltaDel);
357
358 double kappaDel = isothermalCompressibility();
359 double d2rhodp2 = (rhoDel * kappaDel - rho * kappa) / (rhoDel - rho);
360 if (d2rhodp2 > 0.0) {
361 iState = WATER_UNSTABLELIQUID;
362 } else {
363 iState = WATER_UNSTABLEGAS;
364 }
365 delta = deltaSave;
367 }
368 }
369 }
370 return iState;
371}
372
374{
375 double temperature = T_c/tau;
376 double delta_save = delta;
377 // return the critical density if we are above or even just a little below
378 // the critical temperature. We just don't want to worry about the critical
379 // point at this juncture.
380 if (temperature >= T_c - 0.001) {
381 return Rho_c;
382 }
383 double p = psat_est(temperature);
384 double rho_low = 0.0;
385 double rho_high = 1000;
386 double densSatLiq = density_const(p, WATER_LIQUID);
387 double dens_old = densSatLiq;
388 delta = dens_old / Rho_c;
390 double dpdrho_old = dpdrho();
391 if (dpdrho_old > 0.0) {
392 rho_high = std::min(dens_old, rho_high);
393 } else {
394 rho_low = std::max(rho_low, dens_old);
395 }
396 double dens_new = densSatLiq* (1.0001);
397 delta = dens_new / Rho_c;
399 double dpdrho_new = dpdrho();
400 if (dpdrho_new > 0.0) {
401 rho_high = std::min(dens_new, rho_high);
402 } else {
403 rho_low = std::max(rho_low, dens_new);
404 }
405 bool conv = false;
406
407 for (int it = 0; it < 50; it++) {
408 double slope = (dpdrho_new - dpdrho_old)/(dens_new - dens_old);
409 if (slope >= 0.0) {
410 slope = std::max(slope, dpdrho_new *5.0/ dens_new);
411 } else {
412 slope = -dpdrho_new;
413 // shouldn't be here for liquid spinodal
414 }
415 double delta_rho = - dpdrho_new / slope;
416 if (delta_rho > 0.0) {
417 delta_rho = std::min(delta_rho, dens_new * 0.1);
418 } else {
419 delta_rho = std::max(delta_rho, - dens_new * 0.1);
420 }
421 double dens_est = dens_new + delta_rho;
422 if (dens_est < rho_low) {
423 dens_est = 0.5 * (rho_low + dens_new);
424 }
425 if (dens_est > rho_high) {
426 dens_est = 0.5 * (rho_high + dens_new);
427 }
428
429 dens_old = dens_new;
430 dpdrho_old = dpdrho_new;
431 dens_new = dens_est;
432
433 delta = dens_new / Rho_c;
435 dpdrho_new = dpdrho();
436 if (dpdrho_new > 0.0) {
437 rho_high = std::min(dens_new, rho_high);
438 } else if (dpdrho_new < 0.0) {
439 rho_low = std::max(rho_low, dens_new);
440 } else {
441 conv = true;
442 break;
443 }
444
445 if (fabs(dpdrho_new) < 1.0E-5) {
446 conv = true;
447 break;
448 }
449 }
450
451 if (!conv) {
452 throw CanteraError("WaterPropsIAPWS::densSpinodalWater",
453 "convergence failure");
454 }
455 // Restore the original delta
456 delta = delta_save;
458 return dens_new;
459}
460
462{
463 double temperature = T_c/tau;
464 double delta_save = delta;
465 // return the critical density if we are above or even just a little below
466 // the critical temperature. We just don't want to worry about the critical
467 // point at this juncture.
468 if (temperature >= T_c - 0.001) {
469 return Rho_c;
470 }
471 double p = psat_est(temperature);
472 double rho_low = 0.0;
473 double rho_high = 1000;
474 double densSatGas = density_const(p, WATER_GAS);
475 double dens_old = densSatGas;
476 delta = dens_old / Rho_c;
478 double dpdrho_old = dpdrho();
479 if (dpdrho_old < 0.0) {
480 rho_high = std::min(dens_old, rho_high);
481 } else {
482 rho_low = std::max(rho_low, dens_old);
483 }
484 double dens_new = densSatGas * (0.99);
485 delta = dens_new / Rho_c;
487 double dpdrho_new = dpdrho();
488 if (dpdrho_new < 0.0) {
489 rho_high = std::min(dens_new, rho_high);
490 } else {
491 rho_low = std::max(rho_low, dens_new);
492 }
493 bool conv = false;
494 for (int it = 0; it < 50; it++) {
495 double slope = (dpdrho_new - dpdrho_old)/(dens_new - dens_old);
496 if (slope >= 0.0) {
497 slope = dpdrho_new;
498 // shouldn't be here for gas spinodal
499 } else {
500 slope = std::min(slope, dpdrho_new *5.0 / dens_new);
501
502 }
503 double delta_rho = - dpdrho_new / slope;
504 if (delta_rho > 0.0) {
505 delta_rho = std::min(delta_rho, dens_new * 0.1);
506 } else {
507 delta_rho = std::max(delta_rho, - dens_new * 0.1);
508 }
509 double dens_est = dens_new + delta_rho;
510 if (dens_est < rho_low) {
511 dens_est = 0.5 * (rho_low + dens_new);
512 }
513 if (dens_est > rho_high) {
514 dens_est = 0.5 * (rho_high + dens_new);
515 }
516
517 dens_old = dens_new;
518 dpdrho_old = dpdrho_new;
519 dens_new = dens_est;
520 delta = dens_new / Rho_c;
522 dpdrho_new = dpdrho();
523 if (dpdrho_new < 0.0) {
524 rho_high = std::min(dens_new, rho_high);
525 } else if (dpdrho_new > 0.0) {
526 rho_low = std::max(rho_low, dens_new);
527 } else {
528 conv = true;
529 break;
530 }
531
532 if (fabs(dpdrho_new) < 1.0E-5) {
533 conv = true;
534 break;
535 }
536 }
537
538 if (!conv) {
539 throw CanteraError("WaterPropsIAPWS::densSpinodalSteam",
540 "convergence failure");
541 }
542 // Restore the original delta
543 delta = delta_save;
545 return dens_new;
546}
547
548void WaterPropsIAPWS::setState_TD(double temperature, double rho)
549{
550 calcDim(temperature, rho);
552}
553
555{
556 return m_phi.gibbs_RT() * R_water * T_c / tau;
557}
558
560{
561 return m_phi.enthalpy_RT() * R_water * T_c / tau;
562}
563
565{
566 return m_phi.intEnergy_RT() * R_water * T_c / tau;
567}
568
570{
571 return m_phi.entropy_R() * R_water;
572}
573
575{
576 return m_phi.cv_R() * R_water;
577}
578
580{
581 return m_phi.cp_R() * R_water;
582}
583
584}
Headers for a class for calculating the equation of state of water from the IAPWS 1995 Formulation ba...
Base class for exceptions thrown by Cantera classes.
double coeffThermExp() const
Returns the coefficient of thermal expansion.
double densSpinodalSteam() const
Return the value of the density at the water spinodal point (on the gas side) for the current tempera...
double density() const
Returns the density (kg m-3)
void corr(double temperature, double pressure, double &densLiq, double &densGas, double &delGRT)
Utility routine in the calculation of the saturation pressure.
double pressure() const
Calculates the pressure (Pascals), given the current value of the temperature and density.
void corr1(double temperature, double pressure, double &densLiq, double &densGas, double &pcorr)
Utility routine in the calculation of the saturation pressure.
double gibbs_mass() const
Get the Gibbs free energy (J/kg) at the current temperature and density.
double isothermalCompressibility() const
Returns the coefficient of isothermal compressibility for the state of the object.
double psat(double temperature, int waterState=WATER_LIQUID)
This function returns the saturation pressure given the temperature as an input parameter,...
double temperature() const
Returns the temperature (Kelvin)
double cv_mass() const
Get the constant volume heat capacity (J/kg/K) at the current temperature and density.
double entropy_mass() const
Get the entropy (J/kg/K) at the current temperature and density.
double densSpinodalWater() const
Return the value of the density at the water spinodal point (on the liquid side) for the current temp...
double delta
Dimensionless density, delta = rho / rho_c.
double cp_mass() const
Get the constant pressure heat capacity (J/kg/K) at the current temperature and density.
double intEnergy_mass() const
Get the internal energy (J/kg) at the current temperature and density.
double coeffPresExp() const
Returns the isochoric pressure derivative wrt temperature.
int iState
Current state of the system.
void setState_TD(double temperature, double rho)
Set the internal state of the object wrt temperature and density.
double tau
Dimensionless temperature, tau = T_C / T.
WaterPropsIAPWSphi m_phi
pointer to the underlying object that does the calculations.
double density_const(double pressure, int phase=-1, double rhoguess=-1.0) const
Calculates the density given the temperature and the pressure, and a guess at the density,...
void calcDim(double temperature, double rho)
Calculate the dimensionless temp and rho and store internally.
int phaseState(bool checkState=false) const
Returns the Phase State flag for the current state of the object.
double dpdrho() const
Returns the value of dp / drho at constant T for the state of the object.
double psat_est(double temperature) const
This function returns an estimated value for the saturation pressure.
double enthalpy_mass() const
Get the enthalpy (J/kg) at the current temperature and density.
double gibbs_RT() const
Calculate the dimensionless Gibbs free energy.
double cv_R() const
Calculate the dimensionless constant volume heat capacity, Cv/R.
double intEnergy_RT() const
Calculate the dimensionless internal energy, u/RT.
double enthalpy_RT() const
Calculate the dimensionless enthalpy, h/RT.
double entropy_R() const
Calculate the dimensionless entropy, s/R.
double phiR() const
Calculate Equation 6.6 for phiR, the residual part of the dimensionless Helmholtz free energy.
double pressureM_rhoRT(double tau, double delta)
Calculate the dimensionless pressure at tau and delta;.
double dimdpdrho(double tau, double delta)
Dimensionless derivative of p wrt rho at constant T.
double dfind(double p_red, double tau, double deltaGuess)
This function computes the reduced density, given the reduced pressure and the reduced temperature,...
double dimdpdT(double tau, double delta)
Dimensionless derivative of p wrt T at constant rho.
double cp_R() const
Calculate the dimensionless constant pressure heat capacity, Cv/R.
void tdpolycalc(double tau, double delta)
Calculates internal polynomials in tau and delta.
Definitions for the classes that are thrown when Cantera experiences an error condition (also contain...
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
Namespace for the Cantera kernel.
Definition AnyMap.cpp:564
const double Rho_c
Value of the Density at the critical point (kg m-3)
const double T_c
Critical Temperature value (kelvin)
static const double P_c
Critical Pressure (Pascals)
Contains declarations for string manipulation functions within Cantera.