CoolProp  6.6.0
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:
22 
23  public:
25  std::string backend_name(void) {
27  }
28 
29  // REQUIRED BUT NOT USED IN IF97 FUNCTIONS
30  bool using_mole_fractions(void) {
31  return false;
32  };
33  bool using_mass_fractions(void) {
34  return true;
35  }; // But actually it doesn't matter since it is only a pure fluid
36  bool using_volu_fractions(void) {
37  return false;
38  };
39  void set_mole_fractions(const std::vector<CoolPropDbl>& mole_fractions) {
40  throw NotImplementedError("Mole composition has not been implemented.");
41  };
42  void set_mass_fractions(const std::vector<CoolPropDbl>& mass_fractions){}; // Not implemented, but don't throw any errors
43  void set_volu_fractions(const std::vector<CoolPropDbl>& volu_fractions) {
44  throw NotImplementedError("Volume composition has not been implemented.");
45  };
46  const std::vector<CoolPropDbl>& get_mole_fractions(void) {
47  throw NotImplementedError("get_mole_fractions composition has not been implemented.");
48  };
49 
51  bool clear() {
52  // Reset all instances of CachedElement and overwrite
53  // the internal double values with -_HUGE
54  // Default phase condition is no phase imposed
55  // IF97 will make phase/region determination
56  this->_T = -_HUGE;
57  this->_p = -_HUGE;
58  this->_Q = -_HUGE;
59  this->_rhomass.clear();
60  this->_hmass.clear();
61  this->_smass.clear();
62  this->_phase = iphase_not_imposed;
63  return true;
64  };
65 
66  void set_phase() {
67  double epsilon = 3.3e-5; // IAPWS-IF97 RMS saturated pressure inconsistency
68  if ((std::abs(_T - IF97::Tcrit) < epsilon / 10.0) && // RMS temperature inconsistency ~ epsilon/10
69  (std::abs(_p - IF97::Pcrit) < epsilon)) { // within epsilon of [Tcrit,Pcrit]
70  _phase = iphase_critical_point; // at critical point
71  } else if (_T >= IF97::Tcrit) { // to the right of the critical point
72  if (_p >= IF97::Pcrit) { // above the critical point
74  } else { // below the critical point
76  }
77  } else { // to the left of the critical point
78  if (_p >= IF97::Pcrit) { // above the critical point
80  } else { // below critical point
81  double psat = IF97::psat97(_T);
82  if (_p > psat * (1.0 + epsilon)) { // above the saturation curve
84  } else if (_p < psat * (1.0 - epsilon)) { // below the saturation curve
86  } else // exactly on saturation curve (within 1e-4 %)
88  }
89  }
90  };
91 
93 
101  void update(CoolProp::input_pairs input_pair, double value1, double value2) {
102 
103  double H, S, hLmass, hVmass, sLmass, sVmass;
104 
105  clear(); //clear the few cached values we are using
106 
107  switch (input_pair) {
108  case PT_INPUTS:
109  _p = value1;
110  _T = value2;
111  _Q = -1;
112  set_phase();
113  //Two-Phase Check, with PT Inputs:
114  if (_phase == iphase_twophase)
115  throw ValueError(format("Saturation pressure [%g Pa] corresponding to T [%g K] is within 3.3e-3 %% of given p [%Lg Pa]",
116  IF97::psat97(_T), _T, _p));
117  break;
118  case PQ_INPUTS:
119  _p = value1;
120  _Q = value2;
121  if ((_Q < 0) || (_Q > 1)) throw CoolProp::OutOfRangeError("Input vapor quality [Q] must be between 0 and 1");
122  _T = IF97::Tsat97(_p); // ...will throw exception if _P not on saturation curve
124  break;
125  case QT_INPUTS:
126  _Q = value1;
127  _T = value2;
128  if ((_Q < 0) || (_Q > 1)) throw CoolProp::OutOfRangeError("Input vapor quality [Q] must be between 0 and 1");
129  _p = IF97::psat97(_T); // ...will throw exception if _P not on saturation curve
131  break;
132  case HmolarP_INPUTS:
133  // IF97 is mass based so convert hmolar input to hmass
134  _hmass = value1 / molar_mass(); // Convert to mass basis : (J/mol) / (kg/mol) = J/kg
135  _p = value2;
136  // Fall thru to mass basis inputs
137  case HmassP_INPUTS:
138  if (!(_hmass)) _hmass = value1; // don't set if already set above
139  _p = value2;
140  _T = IF97::T_phmass(_p, _hmass);
141  // ...if in the vapor dome (Region 4), calculate Quality...
142  if (IF97::BackwardRegion(_p, _hmass, IF97_HMASS) == 4) {
143  H = _hmass;
144  hVmass = IF97::hvap_p(_p);
145  hLmass = IF97::hliq_p(_p);
146  _Q = (H - hLmass) / (hVmass - hLmass);
148  } else {
149  _Q = -1;
150  set_phase();
151  };
152  break;
153  case PSmolar_INPUTS:
154  // IF97 is mass based so convert smolar input to smass
155  _p = value1;
156  _smass = value2 / molar_mass(); // Convert to mass basis : (J/mol-K) / (kg/mol) = J/kg-K
157  // Fall thru to mass basis inputs
158  case PSmass_INPUTS:
159  _p = value1;
160  if (!(_smass)) _smass = value2;
161  _T = IF97::T_psmass(_p, _smass);
162  if (IF97::BackwardRegion(_p, _smass, IF97_SMASS) == 4) {
163  S = _smass;
164  sVmass = IF97::svap_p(_p);
165  sLmass = IF97::sliq_p(_p);
166  _Q = (S - sLmass) / (sVmass - sLmass);
168  } else {
169  _Q = -1;
170  set_phase();
171  };
172  break;
173  case HmolarSmolar_INPUTS:
174  // IF97 is mass based so convert smolar input to smass
175  _hmass = value1 / molar_mass(); // Convert to mass basis : (J/mol) / (kg/mol) = J/kg
176  _smass = value2 / molar_mass(); // Convert to mass basis : (J/mol-K) / (kg/mol) = J/kg-K
177  // Fall thru to mass basis inputs
178  case HmassSmass_INPUTS:
179  _hmass = value1;
180  _smass = value2;
181  _p = IF97::p_hsmass(_hmass, _smass);
182  _T = IF97::T_phmass(_p, _hmass);
183  // ...if in the vapor dome (Region 4), calculate Quality...
184  if (IF97::BackwardRegion(_p, _hmass, IF97_HMASS) == 4) {
185  H = _hmass;
186  hVmass = IF97::hvap_p(_p);
187  hLmass = IF97::hliq_p(_p);
188  _Q = (H - hLmass) / (hVmass - hLmass);
190  } else {
191  _Q = -1;
192  set_phase();
193  };
194  break;
195  default:
196  throw ValueError("This pair of inputs is not yet supported");
197  }
198  };
199 
200  /* We have to override some of the functions from the AbstractState.
201  * IF97 is only mass-based and does not support conversion
202  * from mass- to molar-specific quantities.
203  */
204  // ************************************************************************* //
205  // Basic Thermodynamic Functions //
206  // ************************************************************************* //
207  //
208  double calc_SatLiquid(parameters iCalc) {
209  switch (iCalc) {
210  case iDmass:
211  return IF97::rholiq_p(_p);
212  break;
213  case iHmass:
214  return IF97::hliq_p(_p);
215  break;
216  case iSmass:
217  return IF97::sliq_p(_p);
218  break;
219  case iCpmass:
220  return IF97::cpliq_p(_p);
221  break;
222  case iCvmass:
223  return IF97::cvliq_p(_p);
224  break;
225  case iUmass:
226  return IF97::uliq_p(_p);
227  break;
228  case ispeed_sound:
229  return IF97::speed_soundliq_p(_p);
230  break;
231  case iviscosity:
232  return IF97::viscliq_p(_p);
233  break;
234  case iconductivity:
235  return IF97::tcondliq_p(_p);
236  break;
237  case isurface_tension:
238  return IF97::sigma97(_T);
239  break;
240  case iPrandtl:
241  return IF97::prandtlliq_p(_p);
242  break;
243  default:
244  return -_HUGE;
245  };
246  }
247  double calc_SatVapor(parameters iCalc) {
248  switch (iCalc) {
249  case iDmass:
250  return IF97::rhovap_p(_p);
251  break;
252  case iHmass:
253  return IF97::hvap_p(_p);
254  break;
255  case iSmass:
256  return IF97::svap_p(_p);
257  break;
258  case iCpmass:
259  return IF97::cpvap_p(_p);
260  break;
261  case iCvmass:
262  return IF97::cvvap_p(_p);
263  break;
264  case iUmass:
265  return IF97::uvap_p(_p);
266  break;
267  case ispeed_sound:
268  return IF97::speed_soundvap_p(_p);
269  break;
270  case iviscosity:
271  return IF97::viscvap_p(_p);
272  break;
273  case iconductivity:
274  return IF97::tcondvap_p(_p);
275  break;
276  case isurface_tension:
277  return IF97::sigma97(_T);
278  break;
279  case iPrandtl:
280  return IF97::prandtlvap_p(_p);
281  break;
282  default:
283  return -_HUGE;
284  };
285  }
286  double calc_Flash(parameters iCalc) {
287  switch (_phase) {
288  case iphase_twophase: // In saturation envelope
289  if (std::abs(_Q) < 1e-10) {
290  return calc_SatLiquid(iCalc); // bubble point (Q == 0) on Sat. Liquid curve
291  } else if (std::abs(_Q - 1) < 1e-10) {
292  return calc_SatVapor(iCalc); // dew point (Q == 1) on Sat. Vapor curve
293  } else { // else "inside" bubble ( 0 < Q < 1 )
294  switch (iCalc) {
295  case iDmass:
296  // Density is an inverse phase weighted property, since it's the inverse of specific volume
297  return 1.0 / (_Q / calc_SatVapor(iDmass) + (1.0 - _Q) / calc_SatLiquid(iDmass));
298  break;
299  case iCpmass:
300  throw NotImplementedError(format("Isobaric Specific Heat not valid in two phase region"));
301  break;
302  case iCvmass:
303  throw NotImplementedError(format("Isochoric Specific Heat not valid in two phase region"));
304  break;
305  case ispeed_sound:
306  throw NotImplementedError(format("Speed of Sound not valid in two phase region"));
307  break;
308  case iviscosity:
309  throw NotImplementedError(format("Viscosity not valid in two phase region"));
310  break;
311  case iconductivity:
312  throw NotImplementedError(format("Viscosity not valid in two phase region"));
313  break;
314  case isurface_tension:
315  return IF97::sigma97(_T);
316  break; // Surface Tension is not a phase-weighted property
317  case iPrandtl:
318  throw NotImplementedError(format("Prandtl number is not valid in two phase region"));
319  break;
320  default:
321  return _Q * calc_SatVapor(iCalc) + (1.0 - _Q) * calc_SatLiquid(iCalc); // Phase weighted combination
322  };
323  }
324  break;
325  default: // Outside saturation envelope (iphase_not_imposed), let IF97 determine phase/region
326  switch (iCalc) {
327  case iDmass:
328  return IF97::rhomass_Tp(_T, _p);
329  break;
330  case iHmass:
331  return IF97::hmass_Tp(_T, _p);
332  break;
333  case iSmass:
334  return IF97::smass_Tp(_T, _p);
335  break;
336  case iCpmass:
337  return IF97::cpmass_Tp(_T, _p);
338  break;
339  case iCvmass:
340  return IF97::cvmass_Tp(_T, _p);
341  break;
342  case iUmass:
343  return IF97::umass_Tp(_T, _p);
344  break;
345  case ispeed_sound:
346  return IF97::speed_sound_Tp(_T, _p);
347  break;
348  case iviscosity:
349  return IF97::visc_Tp(_T, _p);
350  break;
351  case iconductivity:
352  return IF97::tcond_Tp(_T, _p);
353  break;
354  case isurface_tension:
355  throw NotImplementedError(format("Surface Tension is only valid within the two phase region; Try PQ or QT inputs"));
356  break;
357  case iPrandtl:
358  return IF97::prandtl_Tp(_T, _p);
359  break;
360  default:
361  throw NotImplementedError(format("Output variable not implemented in IF97 Backend"));
362  };
363  }
364  }
366  double rhomass(void) {
367  return calc_rhomass();
368  };
369  double calc_rhomass(void) {
370  return calc_Flash(iDmass);
371  };
373  double rhomolar(void) {
374  return calc_rhomolar();
375  };
376  double calc_rhomolar(void) {
377  return rhomass() / molar_mass();
378  };
379 
381  double hmass(void) {
382  return calc_hmass();
383  };
384  double calc_hmass(void) {
385  return calc_Flash(iHmass);
386  };
388  double hmolar(void) {
389  return calc_hmolar();
390  };
391  double calc_hmolar(void) {
392  return hmass() * molar_mass();
393  };
394 
396  double smass(void) {
397  return calc_smass();
398  };
399  double calc_smass(void) {
400  return calc_Flash(iSmass);
401  };
403  double smolar(void) {
404  return calc_smolar();
405  };
406  double calc_smolar(void) {
407  return smass() * molar_mass();
408  };
409 
411  double umass(void) {
412  return calc_umass();
413  };
414  double calc_umass(void) {
415  return calc_Flash(iUmass);
416  };
418  double umolar(void) {
419  return calc_umolar();
420  };
421  double calc_umolar(void) {
422  return umass() * molar_mass();
423  };
424 
426  double cpmass(void) {
427  return calc_cpmass();
428  };
429  double calc_cpmass(void) {
430  return calc_Flash(iCpmass);
431  };
433  double cpmolar(void) {
434  return calc_cpmolar();
435  };
436  double calc_cpmolar(void) {
437  return cpmass() * molar_mass();
438  };
439 
441  double cvmass(void) {
442  return calc_cvmass();
443  };
444  double calc_cvmass(void) {
445  return calc_Flash(iCvmass);
446  };
448  double cvmolar(void) {
449  return calc_cvmolar();
450  };
451  double calc_cvmolar(void) {
452  return cvmass() * molar_mass();
453  };
454 
456  double speed_sound(void) {
457  return calc_speed_sound();
458  };
459  double calc_speed_sound(void) {
460  return calc_Flash(ispeed_sound);
461  };
462 
463  // Return the phase
465  return _phase;
466  };
467 
468  //
469  // ************************************************************************* //
470  // Trivial Functions //
471  // ************************************************************************* //
472  //
474  double calc_Ttriple(void) {
475  return IF97::get_Ttrip();
476  };
478  double calc_p_triple(void) {
479  return IF97::get_ptrip();
480  };
482  double calc_T_critical(void) {
483  return IF97::get_Tcrit();
484  };
486  double calc_p_critical(void) {
487  return IF97::get_pcrit();
488  };
491  double calc_gas_constant(void) {
492  return IF97::get_Rgas() * molar_mass();
493  };
495  double calc_molar_mass(void) {
496  return IF97::get_MW();
497  };
499  double calc_acentric_factor(void) {
500  return IF97::get_Acentric();
501  };
503  // TODO: May want to adjust this based on _T, since Region 5
504  // is limited to 50 MPa, instead of 100 MPa elsewhere.
505  double calc_pmax(void) {
506  return IF97::get_Pmax();
507  };
510  double calc_Tmax(void) {
511  return IF97::get_Tmax();
512  };
514  double calc_Tmin(void) {
515  return IF97::get_Tmin();
516  };
519  double rhomolar_critical(void) {
520  return calc_rhomass_critical() / molar_mass();
521  }
522  double rhomass_critical(void) {
523  return calc_rhomass_critical();
524  }
525  // Overwrite the virtual calc_ functions for density
526  double calc_rhomolar_critical(void) {
527  return rhomass_critical() / molar_mass();
528  };
529  double calc_rhomass_critical(void) {
530  return IF97::get_rhocrit();
531  };
532  //
533  // ************************************************************************* //
534  // Saturation Functions //
535  // ************************************************************************* //
536  //
537  double calc_pressure(void) {
538  return _p;
539  };
540  //
541  // ************************************************************************* //
542  // Transport Property Functions //
543  // ************************************************************************* //
544  //
545  // Return viscosity in [Pa-s]
546  double viscosity(void) {
547  return calc_viscosity();
548  };
549  double calc_viscosity(void) {
550  return calc_Flash(iviscosity);
551  };
552  // Return thermal conductivity in [W/m-K]
553  double conductivity(void) {
554  return calc_conductivity();
555  };
556  double calc_conductivity(void) {
557  return calc_Flash(iconductivity);
558  };
559  // Return surface tension in [N/m]
560  double surface_tension(void) {
561  return calc_surface_tension();
562  };
563  double calc_surface_tension(void) {
565  };
566  // Return Prandtl number (mu*Cp/k) [dimensionless]
567  double Prandtl(void) {
568  return calc_Flash(iPrandtl);
569  };
570 };
571 
572 } /* namespace CoolProp */
573 #endif /* IF97BACKEND_H_ */