CoolProp 6.8.1dev
An open-source fluid property and humid air property database
IF97Backend.h
Go to the documentation of this file.
1
2#ifndef IF97BACKEND_H_
3#define IF97BACKEND_H_
4
5#include "DataStructures.h"
6#include "externals/IF97/IF97.h"
7#include "AbstractState.h"
8#include "Exceptions.h"
9#include <vector>
10#include <cmath>
11
12namespace CoolProp {
13
15{
16
17 protected:
21 _reverse; // Also need a way to flag using the IF97 reverse calcs for h(p,s) and s(p,h)
23
24 public:
26 std::string backend_name(void) {
28 }
29
30 // REQUIRED BUT NOT USED IN IF97 FUNCTIONS
32 return false;
33 };
35 return true;
36 }; // But actually it doesn't matter since it is only a pure fluid
38 return false;
39 };
40 void set_mole_fractions(const std::vector<CoolPropDbl>& mole_fractions) {
41 throw NotImplementedError("Mole composition has not been implemented.");
42 };
43 void set_mass_fractions(const std::vector<CoolPropDbl>& mass_fractions){}; // Not implemented, but don't throw any errors
44 void set_volu_fractions(const std::vector<CoolPropDbl>& volu_fractions) {
45 throw NotImplementedError("Volume composition has not been implemented.");
46 };
47 const std::vector<CoolPropDbl>& get_mole_fractions(void) {
48 throw NotImplementedError("get_mole_fractions composition has not been implemented.");
49 };
50
52 bool clear() {
53 // Reset all instances of CachedElement and overwrite
54 // the internal double values with -_HUGE
55 // Default phase condition is no phase imposed
56 // IF97 will make phase/region determination
57 this->_T = -_HUGE;
58 this->_p = -_HUGE;
59 this->_Q = -_HUGE;
60 this->_rhomass.clear();
61 this->_hmass.clear();
62 this->_smass.clear();
63 this->_reverse.clear();
65 return true;
66 };
67
68 // Set phase based on cached values in _p & _T or, if reverse, from (_p,_smass) or (_p,_hmass)
69 void set_phase() {
70 double epsilon = 3.3e-5; // IAPWS-IF97 RMS saturated pressure inconsistency
71 if (!_reverse) {
72 if ((std::abs(_T - IF97::Tcrit) < epsilon / 10.0) && // RMS temperature inconsistency ~ epsilon/10
73 (std::abs(_p - IF97::Pcrit) < epsilon)) { // within epsilon of [Tcrit,Pcrit]
74 _phase = iphase_critical_point; // at critical point
75 } else if (_T > IF97::Tcrit) { // to the right of the critical point
76 if (_p > IF97::Pcrit) { // above the critical point pressure
78 } else { // below the critical point pressure
80 }
81 } else { // to the left of the critical point
82 if (_p > IF97::Pcrit) { // above the critical point pressure
84 } else { // at or below critical point pressure
85 double psat = IF97::psat97(_T);
86 if (_p > psat * (1.0 + epsilon)) { // above the saturation curve
88 } else if (_p < psat * (1.0 - epsilon)) { // below the saturation curve
90 } else // exactly on saturation curve (within 1e-4 %)
92 }
93 }
94 } else { // Backwards: Determine phase from _p & _smass or _hmass
95 int IF97Region;
96 if (_smass) { // Get IF97 Region
97 IF97Region = IF97::BackwardRegion(_p, _smass, IF97_SMASS); // using p, s
98 } else {
99 IF97Region = IF97::BackwardRegion(_p, _hmass, IF97_HMASS); // using p, h
100 }
101 switch (IF97Region) { // Convert IF97 Region to CP iPhase
102 case 1: // IF97::REGION1
103 if (_p <= IF97::Pcrit) {
105 } else {
107 };
108 break;
109 case 2: // IF97::REGION2
110 if (_T <= IF97::Tcrit) {
112 } else {
114 };
115 break;
116 case 3: // IF97::REGION3
117 if (_T < IF97::Tsat97(_p)) {
118 if (_p <= IF97::Pcrit) {
120 } else {
122 };
123 } else {
124 if (_T <= IF97::Tcrit) {
126 } else {
128 }
129 };
130 break;
131 case 4: // IF97::REGION4 (Saturation
132 _phase = iphase_twophase; // already handled but here for completeness
133 break;
134 case 5: // IF97::REGION5 or O.B.
135 default:
136 throw CoolProp::OutOfRangeError("Outside of IF97 Reverse Function Bounds");
137 break;
138 };
139 }
140 };
141
143
151 void update(CoolProp::input_pairs input_pair, double value1, double value2) {
152
153 double H, S, hLmass, hVmass, sLmass, sVmass;
154
155 clear(); //clear the few cached values we are using
156
157 switch (input_pair) {
158 case PT_INPUTS:
159 _p = value1;
160 _T = value2;
161 _Q = -1;
162 set_phase();
163 //Two-Phase Check, with PT Inputs:
164 if (_phase == iphase_twophase)
165 throw ValueError(format("Saturation pressure [%g Pa] corresponding to T [%g K] is within 3.3e-3 %% of given p [%Lg Pa]",
166 IF97::psat97(_T), _T, _p));
167 break;
168 case PQ_INPUTS:
169 _p = value1;
170 _Q = value2;
171 if ((_Q < 0) || (_Q > 1)) throw CoolProp::OutOfRangeError("Input vapor quality [Q] must be between 0 and 1");
172 _T = IF97::Tsat97(_p); // ...will throw exception if _P not on saturation curve
174 break;
175 case QT_INPUTS:
176 _Q = value1;
177 _T = value2;
178 if ((_Q < 0) || (_Q > 1)) throw CoolProp::OutOfRangeError("Input vapor quality [Q] must be between 0 and 1");
179 _p = IF97::psat97(_T); // ...will throw exception if _P not on saturation curve
181 break;
182 case HmolarP_INPUTS:
183 // IF97 is mass based so convert hmolar input to hmass
184 _hmass = value1 / molar_mass(); // Convert to mass basis : (J/mol) / (kg/mol) = J/kg
185 _p = value2;
186 // Fall thru to mass basis inputs
187 case HmassP_INPUTS:
188 if (!(_hmass)) _hmass = value1; // don't set if already set above
189 _p = value2;
190 _reverse = 1.0; // only ever set for HP and SP Inputs
191 _T = IF97::T_phmass(_p, _hmass);
192 _Q = -1; // Default if not set in 2-Phase Region
193 // ...if in the vapor dome (Region 4), calculate Quality...
194 if (IF97::BackwardRegion(_p, _hmass, IF97_HMASS) == 4) {
195 H = _hmass;
196 hVmass = IF97::hvap_p(_p);
197 hLmass = IF97::hliq_p(_p);
198 _Q = std::min(1.0, std::max(0.0, (H - hLmass) / (hVmass - hLmass))); //bound between 0 and 1
200 } else {
201 set_phase();
202 };
203 break;
204 case PSmolar_INPUTS:
205 // IF97 is mass based so convert smolar input to smass
206 _p = value1;
207 _smass = value2 / molar_mass(); // Convert to mass basis : (J/mol-K) / (kg/mol) = J/kg-K
208 // Fall thru to mass basis inputs
209 case PSmass_INPUTS:
210 _p = value1;
211 if (!(_smass)) _smass = value2;
212 _reverse = 1.0;
213 _T = IF97::T_psmass(_p, _smass);
214 _Q = -1;
215 // ...if in the vapor dome (Region 4), calculate Quality...
216 if (IF97::BackwardRegion(_p, _smass, IF97_SMASS) == 4) {
217 S = _smass;
218 sVmass = IF97::svap_p(_p);
219 sLmass = IF97::sliq_p(_p);
220 _Q = std::min(1.0, std::max(0.0, (S - sLmass) / (sVmass - sLmass))); //bound between 0 and 1
222 } else {
223 set_phase();
224 };
225 break;
227 // IF97 is mass based so convert smolar input to smass
228 _hmass = value1 / molar_mass(); // Convert to mass basis : (J/mol) / (kg/mol) = J/kg
229 _smass = value2 / molar_mass(); // Convert to mass basis : (J/mol-K) / (kg/mol) = J/kg-K
230 // Fall thru to mass basis inputs
232 _hmass = value1;
233 _smass = value2;
234 _p = IF97::p_hsmass(_hmass, _smass);
235 _T = IF97::T_phmass(_p, _hmass);
236 // ...if in the vapor dome (Region 4), calculate Quality...
237 if (IF97::BackwardRegion(_p, _hmass, IF97_HMASS) == 4) {
238 H = _hmass;
239 hVmass = IF97::hvap_p(_p);
240 hLmass = IF97::hliq_p(_p);
241 _Q = std::min(1.0, std::max(0.0, (H - hLmass) / (hVmass - hLmass))); //bount between 0 and 1
243 } else {
244 _Q = -1;
245 set_phase();
246 };
247 break;
248 default:
249 throw ValueError("This pair of inputs is not yet supported");
250 }
251 };
252
253 /* We have to override some of the functions from the AbstractState.
254 * IF97 is only mass-based and does not support conversion
255 * from mass- to molar-specific quantities.
256 */
257 // ************************************************************************* //
258 // Basic Thermodynamic Functions //
259 // ************************************************************************* //
260 //
262 switch (iCalc) {
263 case iDmass:
264 return IF97::rholiq_p(_p);
265 break;
266 case iHmass:
267 return IF97::hliq_p(_p);
268 break;
269 case iSmass:
270 return IF97::sliq_p(_p);
271 break;
272 case iCpmass:
273 return IF97::cpliq_p(_p);
274 break;
275 case iCvmass:
276 return IF97::cvliq_p(_p);
277 break;
278 case iUmass:
279 return IF97::uliq_p(_p);
280 break;
281 case ispeed_sound:
282 return IF97::speed_soundliq_p(_p);
283 break;
284 case iviscosity:
285 return IF97::viscliq_p(_p);
286 break;
287 case iconductivity:
288 return IF97::tcondliq_p(_p);
289 break;
290 case isurface_tension:
291 return IF97::sigma97(_T);
292 break;
293 case iPrandtl:
294 return IF97::prandtlliq_p(_p);
295 break;
296 default:
297 return -_HUGE;
298 };
299 }
300 double calc_SatVapor(parameters iCalc) {
301 switch (iCalc) {
302 case iDmass:
303 return IF97::rhovap_p(_p);
304 break;
305 case iHmass:
306 return IF97::hvap_p(_p);
307 break;
308 case iSmass:
309 return IF97::svap_p(_p);
310 break;
311 case iCpmass:
312 return IF97::cpvap_p(_p);
313 break;
314 case iCvmass:
315 return IF97::cvvap_p(_p);
316 break;
317 case iUmass:
318 return IF97::uvap_p(_p);
319 break;
320 case ispeed_sound:
321 return IF97::speed_soundvap_p(_p);
322 break;
323 case iviscosity:
324 return IF97::viscvap_p(_p);
325 break;
326 case iconductivity:
327 return IF97::tcondvap_p(_p);
328 break;
329 case isurface_tension:
330 return IF97::sigma97(_T);
331 break;
332 case iPrandtl:
333 return IF97::prandtlvap_p(_p);
334 break;
335 default:
336 return -_HUGE;
337 };
338 }
339 double calc_Flash(parameters iCalc) {
340 switch (_phase) {
341 case iphase_twophase: // In saturation envelope
342 if (std::abs(_Q) < 1e-10) {
343 return calc_SatLiquid(iCalc); // bubble point (Q == 0) on Sat. Liquid curve
344 } else if (std::abs(_Q - 1) < 1e-10) {
345 return calc_SatVapor(iCalc); // dew point (Q == 1) on Sat. Vapor curve
346 } else { // else "inside" bubble ( 0 < Q < 1 )
347 switch (iCalc) {
348 case iDmass:
349 // Density is an inverse phase weighted property, since it's the inverse of specific volume
350 return 1.0 / (_Q / calc_SatVapor(iDmass) + (1.0 - _Q) / calc_SatLiquid(iDmass));
351 break;
352 case iCpmass:
353 throw NotImplementedError(format("Isobaric Specific Heat not valid in two phase region"));
354 break;
355 case iCvmass:
356 throw NotImplementedError(format("Isochoric Specific Heat not valid in two phase region"));
357 break;
358 case ispeed_sound:
359 throw NotImplementedError(format("Speed of Sound not valid in two phase region"));
360 break;
361 case iviscosity:
362 throw NotImplementedError(format("Viscosity not valid in two phase region"));
363 break;
364 case iconductivity:
365 throw NotImplementedError(format("Viscosity not valid in two phase region"));
366 break;
367 case isurface_tension:
368 return IF97::sigma97(_T);
369 break; // Surface Tension is not a phase-weighted property
370 case iPrandtl:
371 throw NotImplementedError(format("Prandtl number is not valid in two phase region"));
372 break;
373 default:
374 return _Q * calc_SatVapor(iCalc) + (1.0 - _Q) * calc_SatLiquid(iCalc); // Phase weighted combination
375 };
376 }
377 break;
378 default: // Outside saturation envelope (iphase_not_imposed), let IF97 determine phase/region
379 switch (iCalc) {
380 case iDmass:
381 return IF97::rhomass_Tp(_T, _p);
382 break;
383 case iHmass:
384 return IF97::hmass_Tp(_T, _p);
385 break;
386 case iSmass:
387 return IF97::smass_Tp(_T, _p);
388 break;
389 case iCpmass:
390 return IF97::cpmass_Tp(_T, _p);
391 break;
392 case iCvmass:
393 return IF97::cvmass_Tp(_T, _p);
394 break;
395 case iUmass:
396 return IF97::umass_Tp(_T, _p);
397 break;
398 case ispeed_sound:
399 return IF97::speed_sound_Tp(_T, _p);
400 break;
401 case iviscosity:
402 return IF97::visc_Tp(_T, _p);
403 break;
404 case iconductivity:
405 return IF97::tcond_Tp(_T, _p);
406 break;
407 case isurface_tension:
408 throw NotImplementedError(format("Surface Tension is only valid within the two phase region; Try PQ or QT inputs"));
409 break;
410 case iPrandtl:
411 return IF97::prandtl_Tp(_T, _p);
412 break;
413 default:
414 throw NotImplementedError(format("Output variable not implemented in IF97 Backend"));
415 };
416 }
417 }
419 double rhomass(void) {
420 return calc_rhomass();
421 };
422 double calc_rhomass(void) {
423 return calc_Flash(iDmass);
424 };
426 double rhomolar(void) {
427 return calc_rhomolar();
428 };
429 double calc_rhomolar(void) {
430 return rhomass() / molar_mass();
431 };
432
434 double hmass(void) {
435 return calc_hmass();
436 };
437 double calc_hmass(void) {
438 if (_reverse && _smass)
439 return IF97::hmass_psmass(_p, _smass); // Special IF97 function for handling reverse h(p,s) evaluation
440 else
441 return calc_Flash(iHmass);
442 };
444 double hmolar(void) {
445 return calc_hmolar();
446 };
447 double calc_hmolar(void) {
448 return hmass() * molar_mass();
449 };
450
452 double smass(void) {
453 return calc_smass();
454 };
455 double calc_smass(void) {
456 if (_reverse && _hmass)
457 return IF97::smass_phmass(_p, _hmass); // Special IF97 function for handling reverse s(p,h) evaluation
458 else
459 return calc_Flash(iSmass);
460 };
462 double smolar(void) {
463 return calc_smolar();
464 };
465 double calc_smolar(void) {
466 return smass() * molar_mass();
467 };
468
470 double umass(void) {
471 return calc_umass();
472 };
473 double calc_umass(void) {
474 return calc_Flash(iUmass);
475 };
477 double umolar(void) {
478 return calc_umolar();
479 };
480 double calc_umolar(void) {
481 return umass() * molar_mass();
482 };
483
485 double cpmass(void) {
486 return calc_cpmass();
487 };
488 double calc_cpmass(void) {
489 return calc_Flash(iCpmass);
490 };
492 double cpmolar(void) {
493 return calc_cpmolar();
494 };
495 double calc_cpmolar(void) {
496 return cpmass() * molar_mass();
497 };
498
500 double cvmass(void) {
501 return calc_cvmass();
502 };
503 double calc_cvmass(void) {
504 return calc_Flash(iCvmass);
505 };
507 double cvmolar(void) {
508 return calc_cvmolar();
509 };
510 double calc_cvmolar(void) {
511 return cvmass() * molar_mass();
512 };
513
515 double speed_sound(void) {
516 return calc_speed_sound();
517 };
518 double calc_speed_sound(void) {
519 return calc_Flash(ispeed_sound);
520 };
521
522 // Return the phase
524 return _phase;
525 };
526
527 //
528 // ************************************************************************* //
529 // Trivial Functions //
530 // ************************************************************************* //
531 //
533 double calc_Ttriple(void) {
534 return IF97::get_Ttrip();
535 };
537 double calc_p_triple(void) {
538 return IF97::get_ptrip();
539 };
541 double calc_T_critical(void) {
542 return IF97::get_Tcrit();
543 };
545 double calc_p_critical(void) {
546 return IF97::get_pcrit();
547 };
550 double calc_gas_constant(void) {
551 return IF97::get_Rgas() * molar_mass();
552 };
554 double calc_molar_mass(void) {
555 return IF97::get_MW();
556 };
558 double calc_acentric_factor(void) {
559 return IF97::get_Acentric();
560 };
562 // TODO: May want to adjust this based on _T, since Region 5
563 // is limited to 50 MPa, instead of 100 MPa elsewhere.
564 double calc_pmax(void) {
565 return IF97::get_Pmax();
566 };
569 double calc_Tmax(void) {
570 return IF97::get_Tmax();
571 };
573 double calc_Tmin(void) {
574 return IF97::get_Tmin();
575 };
578 double rhomolar_critical(void) {
580 }
581 double rhomass_critical(void) {
582 return calc_rhomass_critical();
583 }
584 // Overwrite the virtual calc_ functions for density
586 return rhomass_critical() / molar_mass();
587 };
589 return IF97::get_rhocrit();
590 };
591 //
592 // ************************************************************************* //
593 // Saturation Functions //
594 // ************************************************************************* //
595 //
596 double calc_pressure(void) {
597 return _p;
598 };
599 //
600 // ************************************************************************* //
601 // Transport Property Functions //
602 // ************************************************************************* //
603 //
604 // Return viscosity in [Pa-s]
605 double viscosity(void) {
606 return calc_viscosity();
607 };
608 double calc_viscosity(void) {
609 return calc_Flash(iviscosity);
610 };
611 // Return thermal conductivity in [W/m-K]
612 double conductivity(void) {
613 return calc_conductivity();
614 };
615 double calc_conductivity(void) {
617 };
618 // Return surface tension in [N/m]
619 double surface_tension(void) {
620 return calc_surface_tension();
621 };
622 double calc_surface_tension(void) {
624 };
625 // Return Prandtl number (mu*Cp/k) [dimensionless]
626 double Prandtl(void) {
627 return calc_Flash(iPrandtl);
628 };
629};
630
631} /* namespace CoolProp */
632#endif /* IF97BACKEND_H_ */