CoolProp  6.6.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 
12 namespace CoolProp {
13 
14 class IF97Backend : public AbstractState
15 {
16 
17  protected:
20  CachedElement _hmass, _rhomass, _smass,
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) {
27  return get_backend_string(IF97_BACKEND);
28  }
29 
30  // REQUIRED BUT NOT USED IN IF97 FUNCTIONS
31  bool using_mole_fractions(void) {
32  return false;
33  };
34  bool using_mass_fractions(void) {
35  return true;
36  }; // But actually it doesn't matter since it is only a pure fluid
37  bool using_volu_fractions(void) {
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();
64  this->_phase = iphase_not_imposed;
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
77  _phase = iphase_supercritical;
78  } else { // below the critical point pressure
79  _phase = iphase_supercritical_gas;
80  }
81  } else { // to the left of the critical point
82  if (_p > IF97::Pcrit) { // above the critical point pressure
83  _phase = iphase_supercritical_liquid;
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
87  _phase = iphase_liquid;
88  } else if (_p < psat * (1.0 - epsilon)) { // below the saturation curve
89  _phase = iphase_gas;
90  } else // exactly on saturation curve (within 1e-4 %)
91  _phase = iphase_twophase;
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) {
104  _phase = iphase_liquid;
105  } else {
106  _phase = iphase_supercritical_liquid;
107  };
108  break;
109  case 2: // IF97::REGION2
110  if (_T <= IF97::Tcrit) {
111  _phase = iphase_gas;
112  } else {
113  _phase = iphase_supercritical_gas;
114  };
115  break;
116  case 3: // IF97::REGION3
117  if (_T < IF97::Tsat97(_p)) {
118  if (_p <= IF97::Pcrit) {
119  _phase = iphase_liquid;
120  } else {
121  _phase = iphase_supercritical_liquid;
122  };
123  } else {
124  if (_T <= IF97::Tcrit) {
125  _phase = iphase_gas;
126  } else {
127  _phase = iphase_supercritical_gas;
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
173  _phase = iphase_twophase;
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
180  _phase = iphase_twophase;
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
199  _phase = iphase_twophase;
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
221  _phase = iphase_twophase;
222  } else {
223  set_phase();
224  };
225  break;
226  case HmolarSmolar_INPUTS:
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
231  case HmassSmass_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
242  _phase = iphase_twophase;
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  //
261  double calc_SatLiquid(parameters iCalc) {
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
523  phases calc_phase(void) {
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) {
579  return calc_rhomass_critical() / molar_mass();
580  }
581  double rhomass_critical(void) {
582  return calc_rhomass_critical();
583  }
584  // Overwrite the virtual calc_ functions for density
585  double calc_rhomolar_critical(void) {
586  return rhomass_critical() / molar_mass();
587  };
588  double calc_rhomass_critical(void) {
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) {
616  return calc_Flash(iconductivity);
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) {
623  return calc_Flash(isurface_tension);
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_ */
Generated on Fri Oct 4 2024 09:06:33 for CoolProp by doxygen 1.9.1