CoolProp  6.6.0
An open-source fluid property and humid air property database
Helmholtz.h
Go to the documentation of this file.
1 
2 #ifndef HELMHOLTZ_H
3 #define HELMHOLTZ_H
4 
5 #include <vector>
6 #include "rapidjson_include.h"
7 //#include "Eigen/Core"
8 #include "time.h"
9 #include "CachedElement.h"
10 #include "Backends/Cubics/GeneralizedCubic.h"
11 #include "crossplatform_shared_ptr.h"
12 
13 namespace CoolProp {
14 
15 // #############################################################################
16 // #############################################################################
17 // #############################################################################
18 // RESIDUAL TERMS
19 // #############################################################################
20 // #############################################################################
21 // #############################################################################
22 
23 #define LIST_OF_DERIVATIVE_VARIABLES \
24  X(alphar) \
25  X(dalphar_ddelta) \
26  X(dalphar_dtau) \
27  X(d2alphar_ddelta2) \
28  X(d2alphar_dtau2) \
29  X(d2alphar_ddelta_dtau) \
30  X(d3alphar_ddelta3) \
31  X(d3alphar_ddelta_dtau2) \
32  X(d3alphar_ddelta2_dtau) \
33  X(d3alphar_dtau3) \
34  X(d4alphar_ddelta4) \
35  X(d4alphar_ddelta3_dtau) \
36  X(d4alphar_ddelta2_dtau2) \
37  X(d4alphar_ddelta_dtau3) \
38  X(d4alphar_dtau4) \
39  X(delta_x_dalphar_ddelta) \
40  X(tau_x_dalphar_dtau) \
41  X(delta2_x_d2alphar_ddelta2) \
42  X(deltatau_x_d2alphar_ddelta_dtau) \
43  X(tau2_x_d2alphar_dtau2)
44 
45 struct HelmholtzDerivatives
46 {
47 #define X(name) CoolPropDbl name;
48  LIST_OF_DERIVATIVE_VARIABLES
49 #undef X
50  CoolPropDbl tau, delta, T_red, rhomolar_red;
51 
52  void reset(CoolPropDbl v) {
53 #define X(name) name = v;
54  LIST_OF_DERIVATIVE_VARIABLES
55 #undef X
56  }
57  HelmholtzDerivatives operator+(const HelmholtzDerivatives& other) const {
58  HelmholtzDerivatives _new;
59 #define X(name) _new.name = name + other.name;
60  LIST_OF_DERIVATIVE_VARIABLES
61 #undef X
62  return _new;
63  }
64  HelmholtzDerivatives operator*(const CoolPropDbl& other) const {
65  HelmholtzDerivatives _new;
66 #define X(name) _new.name = name * other;
67  LIST_OF_DERIVATIVE_VARIABLES
68 #undef X
69  return _new;
70  }
71  HelmholtzDerivatives() {
72  reset(0.0);
73  T_red = _HUGE;
74  rhomolar_red = _HUGE;
75  };
77  double get(std::size_t itau, std::size_t idelta) {
78  if (itau == 0) {
79  if (idelta == 0) {
80  return alphar;
81  } else if (idelta == 1) {
82  return dalphar_ddelta;
83  } else if (idelta == 2) {
84  return d2alphar_ddelta2;
85  } else if (idelta == 3) {
86  return d3alphar_ddelta3;
87  } else if (idelta == 4) {
88  return d4alphar_ddelta4;
89  } else {
90  throw ValueError();
91  }
92  } else if (itau == 1) {
93  if (idelta == 0) {
94  return dalphar_dtau;
95  } else if (idelta == 1) {
96  return d2alphar_ddelta_dtau;
97  } else if (idelta == 2) {
98  return d3alphar_ddelta2_dtau;
99  } else if (idelta == 3) {
100  return d4alphar_ddelta3_dtau;
101  } else {
102  throw ValueError();
103  }
104  } else if (itau == 2) {
105  if (idelta == 0) {
106  return d2alphar_dtau2;
107  } else if (idelta == 1) {
108  return d3alphar_ddelta_dtau2;
109  } else if (idelta == 2) {
110  return d4alphar_ddelta2_dtau2;
111  } else {
112  throw ValueError();
113  }
114  } else if (itau == 3) {
115  if (idelta == 0) {
116  return d3alphar_dtau3;
117  } else if (idelta == 1) {
118  return d4alphar_ddelta_dtau3;
119  } else {
120  throw ValueError();
121  }
122  } else if (itau == 4) {
123  if (idelta == 0) {
124  return d4alphar_dtau4;
125  } else {
126  throw ValueError();
127  }
128  } else {
129  throw ValueError();
130  }
131  }
132 };
133 #undef LIST_OF_DERIVATIVE_VARIABLES
134 
136 
161 class BaseHelmholtzTerm
162 {
163  public:
164  BaseHelmholtzTerm(){};
165  virtual ~BaseHelmholtzTerm(){};
166 
168 
171  virtual CoolPropDbl base(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
172  HelmholtzDerivatives deriv;
173  all(tau, delta, deriv);
174  return deriv.alphar;
175  };
177 
180  virtual CoolPropDbl dTau(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
181  HelmholtzDerivatives deriv;
182  all(tau, delta, deriv);
183  return deriv.dalphar_dtau;
184  };
186 
189  virtual CoolPropDbl dTau2(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
190  HelmholtzDerivatives deriv;
191  all(tau, delta, deriv);
192  return deriv.d2alphar_dtau2;
193  };
195 
198  virtual CoolPropDbl dDelta_dTau(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
199  HelmholtzDerivatives deriv;
200  all(tau, delta, deriv);
201  return deriv.d2alphar_ddelta_dtau;
202  };
204 
207  virtual CoolPropDbl dDelta(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
208  HelmholtzDerivatives deriv;
209  all(tau, delta, deriv);
210  return deriv.dalphar_ddelta;
211  };
213 
216  virtual CoolPropDbl dDelta2(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
217  HelmholtzDerivatives deriv;
218  all(tau, delta, deriv);
219  return deriv.d2alphar_ddelta2;
220  };
222 
225  virtual CoolPropDbl dDelta2_dTau(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
226  HelmholtzDerivatives deriv;
227  all(tau, delta, deriv);
228  return deriv.d3alphar_ddelta2_dtau;
229  };
231 
234  virtual CoolPropDbl dDelta_dTau2(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
235  HelmholtzDerivatives deriv;
236  all(tau, delta, deriv);
237  return deriv.d3alphar_ddelta_dtau2;
238  };
240 
243  virtual CoolPropDbl dTau3(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
244  HelmholtzDerivatives deriv;
245  all(tau, delta, deriv);
246  return deriv.d3alphar_dtau3;
247  };
249 
252  virtual CoolPropDbl dDelta3(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
253  HelmholtzDerivatives deriv;
254  all(tau, delta, deriv);
255  return deriv.d3alphar_ddelta3;
256  };
258 
261  virtual CoolPropDbl dTau4(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
262  HelmholtzDerivatives deriv;
263  all(tau, delta, deriv);
264  return deriv.d4alphar_dtau4;
265  };
266  virtual CoolPropDbl dDelta_dTau3(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
267  HelmholtzDerivatives deriv;
268  all(tau, delta, deriv);
269  return deriv.d4alphar_ddelta_dtau3;
270  };
271  virtual CoolPropDbl dDelta2_dTau2(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
272  HelmholtzDerivatives deriv;
273  all(tau, delta, deriv);
274  return deriv.d4alphar_ddelta2_dtau2;
275  };
276  virtual CoolPropDbl dDelta3_dTau(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
277  HelmholtzDerivatives deriv;
278  all(tau, delta, deriv);
279  return deriv.d4alphar_ddelta3_dtau;
280  };
281  virtual CoolPropDbl dDelta4(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
282  HelmholtzDerivatives deriv;
283  all(tau, delta, deriv);
284  return deriv.d4alphar_ddelta4;
285  };
286 
287  virtual void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw() = 0;
288 };
289 
290 struct ResidualHelmholtzGeneralizedExponentialElement
291 {
293  CoolPropDbl n, d, t;
296  CoolPropDbl c, l_double, omega, m_double, eta1, epsilon1, eta2, epsilon2, beta1, gamma1, beta2, gamma2;
298  int l_int, m_int;
300  bool l_is_int, m_is_int;
301 
302  ResidualHelmholtzGeneralizedExponentialElement() {
303  n = 0;
304  d = 0;
305  t = 0;
306  c = 0;
307  l_double = 0;
308  omega = 0;
309  m_double = 0;
310  eta1 = 0;
311  epsilon1 = 0;
312  eta2 = 0;
313  epsilon2 = 0;
314  beta1 = 0;
315  gamma1 = 0;
316  beta2 = 0;
317  gamma2 = 0;
318  l_int = 0;
319  m_int = 0;
320  l_is_int = false;
321  m_is_int = true;
322  };
323 };
332 class ResidualHelmholtzGeneralizedExponential : public BaseHelmholtzTerm
333 {
334 
335  public:
336  bool delta_li_in_u, tau_mi_in_u, eta1_in_u, eta2_in_u, beta1_in_u, beta2_in_u, finished;
337  std::vector<CoolPropDbl> s;
338  std::size_t N;
339 
340  // These variables are for the exp(u) part
341  // u is given by -c*delta^l_i-omega*tau^m_i-eta1*(delta-epsilon1)-eta2*(delta-epsilon2)^2-beta1*(tau-gamma1)-beta2*(tau-gamma2)^2
342  std::vector<double> n, d, t, c, l_double, omega, m_double, eta1, epsilon1, eta2, epsilon2, beta1, gamma1, beta2, gamma2;
343  // If l_i or m_i are integers, we will store them as integers in order to call pow(double, int) rather than pow(double, double)
344  std::vector<int> l_int, m_int;
345 
346  //Eigen::ArrayXd uE, du_ddeltaE, du_dtauE, d2u_ddelta2E, d2u_dtau2E, d3u_ddelta3E, d3u_dtau3E;
347 
348  std::vector<ResidualHelmholtzGeneralizedExponentialElement> elements;
349  // Default Constructor
350  ResidualHelmholtzGeneralizedExponential()
351  : delta_li_in_u(false), tau_mi_in_u(false), eta1_in_u(false), eta2_in_u(false), beta1_in_u(false), beta2_in_u(false), finished(false), N(0){};
357  void add_Power(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& d, const std::vector<CoolPropDbl>& t,
358  const std::vector<CoolPropDbl>& l) {
359  for (std::size_t i = 0; i < n.size(); ++i) {
360  ResidualHelmholtzGeneralizedExponentialElement el;
361  el.n = n[i];
362  el.d = d[i];
363  el.t = t[i];
364  el.l_double = l[i];
365  el.l_int = (int)el.l_double;
366  if (el.l_double > 0)
367  el.c = 1.0;
368  else
369  el.c = 0.0;
370  elements.push_back(el);
371  }
372  delta_li_in_u = true;
373  };
379  void add_Exponential(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& d, const std::vector<CoolPropDbl>& t,
380  const std::vector<CoolPropDbl>& g, const std::vector<CoolPropDbl>& l) {
381  for (std::size_t i = 0; i < n.size(); ++i) {
382  ResidualHelmholtzGeneralizedExponentialElement el;
383  el.n = n[i];
384  el.d = d[i];
385  el.t = t[i];
386  el.c = g[i];
387  el.l_double = l[i];
388  el.l_int = (int)el.l_double;
389  elements.push_back(el);
390  }
391  delta_li_in_u = true;
392  }
398  void add_Gaussian(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& d, const std::vector<CoolPropDbl>& t,
399  const std::vector<CoolPropDbl>& eta, const std::vector<CoolPropDbl>& epsilon, const std::vector<CoolPropDbl>& beta,
400  const std::vector<CoolPropDbl>& gamma) {
401  for (std::size_t i = 0; i < n.size(); ++i) {
402  ResidualHelmholtzGeneralizedExponentialElement el;
403  el.n = n[i];
404  el.d = d[i];
405  el.t = t[i];
406  el.eta2 = eta[i];
407  el.epsilon2 = epsilon[i];
408  el.beta2 = beta[i];
409  el.gamma2 = gamma[i];
410  elements.push_back(el);
411  }
412  eta2_in_u = true;
413  beta2_in_u = true;
414  };
420  void add_GERG2008Gaussian(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& d, const std::vector<CoolPropDbl>& t,
421  const std::vector<CoolPropDbl>& eta, const std::vector<CoolPropDbl>& epsilon, const std::vector<CoolPropDbl>& beta,
422  const std::vector<CoolPropDbl>& gamma) {
423  for (std::size_t i = 0; i < n.size(); ++i) {
424  ResidualHelmholtzGeneralizedExponentialElement el;
425  el.n = n[i];
426  el.d = d[i];
427  el.t = t[i];
428  el.eta2 = eta[i];
429  el.epsilon2 = epsilon[i];
430  el.eta1 = beta[i];
431  el.epsilon1 = gamma[i];
432  elements.push_back(el);
433  }
434  eta2_in_u = true;
435  eta1_in_u = true;
436  };
442  void add_Lemmon2005(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& d, const std::vector<CoolPropDbl>& t,
443  const std::vector<CoolPropDbl>& l, const std::vector<CoolPropDbl>& m) {
444  for (std::size_t i = 0; i < n.size(); ++i) {
445  ResidualHelmholtzGeneralizedExponentialElement el;
446  el.n = n[i];
447  el.d = d[i];
448  el.t = t[i];
449  el.c = 1.0;
450  el.omega = 1.0;
451  el.l_double = l[i];
452  el.m_double = m[i];
453  el.l_int = (int)el.l_double;
454  el.m_int = (int)el.m_double;
455  elements.push_back(el);
456  }
457  delta_li_in_u = true;
458  tau_mi_in_u = true;
459  };
460 
461  void finish() {
462  n.resize(elements.size());
463  d.resize(elements.size());
464  t.resize(elements.size());
465  c.resize(elements.size());
466  omega.resize(elements.size());
467  l_double.resize(elements.size());
468  l_int.resize(elements.size());
469  m_double.resize(elements.size());
470  m_int.resize(elements.size());
471  epsilon2.resize(elements.size());
472  eta2.resize(elements.size());
473  gamma2.resize(elements.size());
474  beta2.resize(elements.size());
475 
476  for (std::size_t i = 0; i < elements.size(); ++i) {
477  n[i] = elements[i].n;
478  d[i] = elements[i].d;
479  t[i] = elements[i].t;
480  c[i] = elements[i].c;
481  omega[i] = elements[i].omega;
482  l_double[i] = elements[i].l_double;
483  l_int[i] = elements[i].l_int;
484  m_double[i] = elements[i].m_double;
485  m_int[i] = elements[i].m_int;
486  epsilon2[i] = elements[i].epsilon2;
487  eta2[i] = elements[i].eta2;
488  gamma2[i] = elements[i].gamma2;
489  beta2[i] = elements[i].beta2;
490 
491  // See if l is an integer, and store a flag if it is
492  elements[i].l_is_int = (std::abs(static_cast<long>(elements[i].l_double) - elements[i].l_double) < 1e-14);
493  }
494  // uE.resize(elements.size());
495  // du_ddeltaE.resize(elements.size());
496  // du_dtauE.resize(elements.size());
497  // d2u_ddelta2E.resize(elements.size());
498  // d2u_dtau2E.resize(elements.size());
499  // d3u_ddelta3E.resize(elements.size());
500  // d3u_dtau3E.resize(elements.size());
501 
502  finished = true;
503  };
504 
505  void to_json(rapidjson::Value& el, rapidjson::Document& doc);
506 
507  void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
508  //void allEigen(const CoolPropDbl &tau, const CoolPropDbl &delta, HelmholtzDerivatives &derivs) throw();
509 };
510 
511 struct ResidualHelmholtzNonAnalyticElement
512 {
513  CoolPropDbl n, a, b, beta, A, B, C, D;
514 };
515 class ResidualHelmholtzNonAnalytic : public BaseHelmholtzTerm
516 {
517 
518  public:
519  std::size_t N;
520  std::vector<CoolPropDbl> s;
521  std::vector<ResidualHelmholtzNonAnalyticElement> elements;
523  ResidualHelmholtzNonAnalytic() {
524  N = 0;
525  };
527  ~ResidualHelmholtzNonAnalytic(){};
529  ResidualHelmholtzNonAnalytic(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& a, const std::vector<CoolPropDbl>& b,
530  const std::vector<CoolPropDbl>& beta, const std::vector<CoolPropDbl>& A, const std::vector<CoolPropDbl>& B,
531  const std::vector<CoolPropDbl>& C, const std::vector<CoolPropDbl>& D) {
532  N = n.size();
533  s.resize(N);
534  for (std::size_t i = 0; i < n.size(); ++i) {
535  ResidualHelmholtzNonAnalyticElement el;
536  el.n = n[i];
537  el.a = a[i];
538  el.b = b[i];
539  el.beta = beta[i];
540  el.A = A[i];
541  el.B = B[i];
542  el.C = C[i];
543  el.D = D[i];
544  elements.push_back(el);
545  }
546  };
547  void to_json(rapidjson::Value& el, rapidjson::Document& doc);
548  void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
549 };
550 
551 class ResidualHelmholtzGeneralizedCubic : public BaseHelmholtzTerm
552 {
553  protected:
554  shared_ptr<AbstractCubic> m_abstractcubic;
555  std::vector<double> z;
556  public:
557  bool enabled;
558 
560  ResidualHelmholtzGeneralizedCubic() {
561  enabled = false;
562  };
564  ResidualHelmholtzGeneralizedCubic(shared_ptr<AbstractCubic>& ac) : m_abstractcubic(ac) {
565  enabled = true;
566  z = std::vector<double>(1, 1); // Init the vector to [1.0]
567  };
568 
569  void to_json(rapidjson::Value& el, rapidjson::Document& doc);
570  void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
571 };
572 
573 class ResidualHelmholtzGaoB : public BaseHelmholtzTerm
574 {
575  protected:
576  std::vector<double> n, t, d, eta, beta, gamma, epsilon, b;
577 
578  public:
579  bool enabled;
580 
582  ResidualHelmholtzGaoB() {
583  enabled = false;
584  };
585 
587  ResidualHelmholtzGaoB(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& t, const std::vector<CoolPropDbl>& d,
588  const std::vector<CoolPropDbl>& eta, const std::vector<CoolPropDbl>& beta, const std::vector<CoolPropDbl>& gamma,
589  const std::vector<CoolPropDbl>& epsilon, const std::vector<CoolPropDbl>& b)
590  : n(n), t(t), d(d), eta(eta), beta(beta), gamma(gamma), epsilon(epsilon), b(b) {
591  enabled = true;
592  };
593 
594  void to_json(rapidjson::Value& el, rapidjson::Document& doc);
595  void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
596 };
597 
599 class ResidualHelmholtzXiangDeiters : public BaseHelmholtzTerm
600 {
601 
602  public:
603  bool enabled;
604  ResidualHelmholtzGeneralizedExponential phi0, phi1, phi2;
605  CoolPropDbl Tc, pc, rhomolarc, acentric, R, theta;
607  ResidualHelmholtzXiangDeiters() : Tc(_HUGE), pc(_HUGE), rhomolarc(_HUGE), acentric(_HUGE), R(_HUGE), theta(_HUGE) {
608  enabled = false;
609  };
611  ResidualHelmholtzXiangDeiters(const CoolPropDbl Tc, const CoolPropDbl pc, const CoolPropDbl rhomolarc, const CoolPropDbl acentric,
612  const CoolPropDbl R);
613  void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
614 };
615 
616 class ResidualHelmholtzSAFTAssociating : public BaseHelmholtzTerm
617 {
618 
619  protected:
620  double a, m, epsilonbar, vbarn, kappabar;
621 
622  CoolPropDbl Deltabar(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
623  CoolPropDbl dDeltabar_ddelta__consttau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
624  CoolPropDbl d2Deltabar_ddelta2__consttau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
625  CoolPropDbl dDeltabar_dtau__constdelta(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
626  CoolPropDbl d2Deltabar_dtau2__constdelta(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
627  CoolPropDbl d2Deltabar_ddelta_dtau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
628  CoolPropDbl d3Deltabar_dtau3__constdelta(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
629  CoolPropDbl d3Deltabar_ddelta_dtau2(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
630  CoolPropDbl d3Deltabar_ddelta3__consttau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
631  CoolPropDbl d3Deltabar_ddelta2_dtau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
632 
633  CoolPropDbl X(const CoolPropDbl& delta, const CoolPropDbl& Deltabar) const;
634  CoolPropDbl dX_dDeltabar__constdelta(const CoolPropDbl& delta, const CoolPropDbl& Deltabar) const;
635  CoolPropDbl dX_ddelta__constDeltabar(const CoolPropDbl& delta, const CoolPropDbl& Deltabar) const;
636  CoolPropDbl dX_dtau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
637  CoolPropDbl dX_ddelta(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
638  CoolPropDbl d2X_dtau2(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
639  CoolPropDbl d2X_ddeltadtau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
640  CoolPropDbl d2X_ddelta2(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
641 
642  CoolPropDbl d3X_dtau3(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
643  CoolPropDbl d3X_ddelta3(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
644  CoolPropDbl d3X_ddeltadtau2(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
645  CoolPropDbl d3X_ddelta2dtau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
646 
647  CoolPropDbl g(const CoolPropDbl& eta) const;
648  CoolPropDbl dg_deta(const CoolPropDbl& eta) const;
649  CoolPropDbl d2g_deta2(const CoolPropDbl& eta) const;
650  CoolPropDbl d3g_deta3(const CoolPropDbl& eta) const;
651  CoolPropDbl eta(const CoolPropDbl& delta) const;
652 
653  public:
655  ResidualHelmholtzSAFTAssociating() : a(_HUGE), m(_HUGE), epsilonbar(_HUGE), vbarn(_HUGE), kappabar(_HUGE) {
656  disabled = true;
657  };
658 
659  // Constructor
660  ResidualHelmholtzSAFTAssociating(double a, double m, double epsilonbar, double vbarn, double kappabar)
661  : a(a), m(m), epsilonbar(epsilonbar), vbarn(vbarn), kappabar(kappabar) {
662  disabled = false;
663  };
664 
665  bool disabled;
666 
667  //Destructor. No Implementation
668  ~ResidualHelmholtzSAFTAssociating(){};
669 
670  void to_json(rapidjson::Value& el, rapidjson::Document& doc);
671 
672  CoolPropDbl dTau4(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
673  return 1e99;
674  };
675  CoolPropDbl dDelta_dTau3(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
676  return 1e99;
677  };
678  CoolPropDbl dDelta2_dTau2(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
679  return 1e99;
680  };
681  CoolPropDbl dDelta3_dTau(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
682  return 1e99;
683  };
684  CoolPropDbl dDelta4(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
685  return 1e99;
686  };
687 
688  void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& deriv) throw();
689 };
690 
691 class BaseHelmholtzContainer
692 {
693  protected:
694  CachedElement _base, _dDelta, _dTau, _dDelta2, _dTau2, _dDelta_dTau, _dDelta3, _dDelta2_dTau, _dDelta_dTau2, _dTau3;
695  CachedElement _dDelta4, _dDelta3_dTau, _dDelta2_dTau2, _dDelta_dTau3, _dTau4;
696 
697  public:
698  void clear() {
699  _base.clear();
700  _dDelta.clear();
701  _dTau.clear();
702  _dDelta2.clear();
703  _dTau2.clear();
704  _dDelta_dTau.clear();
705  _dDelta3.clear();
706  _dTau3.clear();
707  _dDelta2_dTau.clear();
708  _dDelta_dTau2.clear();
709  _dDelta4.clear();
710  _dDelta3_dTau.clear();
711  _dDelta2_dTau2.clear();
712  _dDelta_dTau3.clear();
713  _dTau4.clear();
714  };
715 
716  virtual void empty_the_EOS() = 0;
717  virtual HelmholtzDerivatives all(const CoolPropDbl tau, const CoolPropDbl delta, bool cache_values) = 0;
718 
719  CoolPropDbl base(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
720  if (!_base || dont_use_cache)
721  return all(tau, delta, false).alphar;
722  else
723  return _base;
724  };
725  CoolPropDbl dDelta(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
726  if (!_dDelta || dont_use_cache)
727  return all(tau, delta, false).dalphar_ddelta;
728  else
729  return _dDelta;
730  };
731  CoolPropDbl dTau(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
732  if (!_dTau || dont_use_cache)
733  return all(tau, delta, false).dalphar_dtau;
734  else
735  return _dTau;
736  };
737  CoolPropDbl dDelta2(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
738  if (!_dDelta2 || dont_use_cache)
739  return all(tau, delta, false).d2alphar_ddelta2;
740  else
741  return _dDelta2;
742  };
743  CoolPropDbl dDelta_dTau(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
744  if (!_dDelta_dTau || dont_use_cache)
745  return all(tau, delta, false).d2alphar_ddelta_dtau;
746  else
747  return _dDelta_dTau;
748  };
749  CoolPropDbl dTau2(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
750  if (!_dTau2 || dont_use_cache)
751  return all(tau, delta, false).d2alphar_dtau2;
752  else
753  return _dTau2;
754  };
755  CoolPropDbl dDelta3(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
756  if (!_dDelta3 || dont_use_cache)
757  return all(tau, delta, false).d3alphar_ddelta3;
758  else
759  return _dDelta3;
760  };
761  CoolPropDbl dDelta2_dTau(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
762  if (!_dDelta2_dTau || dont_use_cache)
763  return all(tau, delta, false).d3alphar_ddelta2_dtau;
764  else
765  return _dDelta2_dTau;
766  };
767  CoolPropDbl dDelta_dTau2(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
768  if (!_dDelta_dTau2 || dont_use_cache)
769  return all(tau, delta, false).d3alphar_ddelta_dtau2;
770  else
771  return _dDelta_dTau2;
772  };
773  CoolPropDbl dTau3(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
774  if (!_dTau3 || dont_use_cache)
775  return all(tau, delta, false).d3alphar_dtau3;
776  else
777  return _dTau3;
778  };
779  CoolPropDbl dDelta4(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
780  return all(tau, delta, false).d4alphar_ddelta4;
781  };
782  CoolPropDbl dDelta3_dTau(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
783  return all(tau, delta, false).d4alphar_ddelta3_dtau;
784  };
785  CoolPropDbl dDelta2_dTau2(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
786  return all(tau, delta, false).d4alphar_ddelta2_dtau2;
787  };
788  CoolPropDbl dDelta_dTau3(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
789  return all(tau, delta, false).d4alphar_ddelta_dtau3;
790  };
791  CoolPropDbl dTau4(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
792  return all(tau, delta, false).d4alphar_dtau4;
793  };
794 };
795 
796 class ResidualHelmholtzContainer : public BaseHelmholtzContainer
797 {
798  public:
799  ResidualHelmholtzNonAnalytic NonAnalytic;
800  ResidualHelmholtzSAFTAssociating SAFT;
801  ResidualHelmholtzGeneralizedExponential GenExp;
802  ResidualHelmholtzGeneralizedCubic cubic;
803  ResidualHelmholtzXiangDeiters XiangDeiters;
804  ResidualHelmholtzGaoB GaoB;
805 
806  void empty_the_EOS() {
807  NonAnalytic = ResidualHelmholtzNonAnalytic();
808  SAFT = ResidualHelmholtzSAFTAssociating();
809  GenExp = ResidualHelmholtzGeneralizedExponential();
810  cubic = ResidualHelmholtzGeneralizedCubic();
811  XiangDeiters = ResidualHelmholtzXiangDeiters();
812  GaoB = ResidualHelmholtzGaoB();
813  };
814 
815  HelmholtzDerivatives all(const CoolPropDbl tau, const CoolPropDbl delta, bool cache_values = false) {
816  HelmholtzDerivatives derivs; // zeros out the elements
817  GenExp.all(tau, delta, derivs);
818  NonAnalytic.all(tau, delta, derivs);
819  SAFT.all(tau, delta, derivs);
820  cubic.all(tau, delta, derivs);
821  XiangDeiters.all(tau, delta, derivs);
822  GaoB.all(tau, delta, derivs);
823  if (cache_values) {
824  _base = derivs.alphar;
825  _dDelta = derivs.dalphar_ddelta;
826  _dTau = derivs.dalphar_dtau;
827  _dDelta2 = derivs.d2alphar_ddelta2;
828  _dTau2 = derivs.d2alphar_dtau2;
829  _dDelta_dTau = derivs.d2alphar_ddelta_dtau;
830  _dDelta3 = derivs.d3alphar_ddelta3;
831  _dTau3 = derivs.d3alphar_dtau3;
832  _dDelta2_dTau = derivs.d3alphar_ddelta2_dtau;
833  _dDelta_dTau2 = derivs.d3alphar_ddelta_dtau2;
834  }
835  return derivs;
836  };
837 };
838 
839 // #############################################################################
840 // #############################################################################
841 // #############################################################################
842 // IDEAL GAS TERMS
843 // #############################################################################
844 // #############################################################################
845 // #############################################################################
846 
848 
853 class IdealHelmholtzLead : public BaseHelmholtzTerm
854 {
855 
856  private:
857  CoolPropDbl a1, a2;
858  bool enabled;
859 
860  public:
861  // Default constructor
862  IdealHelmholtzLead() : a1(_HUGE), a2(_HUGE), enabled(false) {}
863 
864  // Constructor
865  IdealHelmholtzLead(CoolPropDbl a1, CoolPropDbl a2) : a1(a1), a2(a2), enabled(true) {}
866 
867  bool is_enabled() const {
868  return enabled;
869  }
870 
871  void to_json(rapidjson::Value& el, rapidjson::Document& doc) {
872  el.AddMember("type", "IdealHelmholtzLead", doc.GetAllocator());
873  el.AddMember("a1", static_cast<double>(a1), doc.GetAllocator());
874  el.AddMember("a2", static_cast<double>(a2), doc.GetAllocator());
875  };
876 
877  void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
878 };
879 
881 
886 class IdealHelmholtzEnthalpyEntropyOffset : public BaseHelmholtzTerm
887 {
888  private:
889  CoolPropDbl a1, a2; // Use these variables internally
890  std::string reference;
891  bool enabled;
892 
893  public:
894  IdealHelmholtzEnthalpyEntropyOffset() : a1(_HUGE), a2(_HUGE), enabled(false) {}
895 
896  // Constructor
897  IdealHelmholtzEnthalpyEntropyOffset(CoolPropDbl a1, CoolPropDbl a2, const std::string& ref) : a1(a1), a2(a2), reference(ref), enabled(true) {}
898 
899  // Set the values in the class
900  void set(CoolPropDbl a1, CoolPropDbl a2, const std::string& ref) {
901  // If it doesn't already exist, just set the values
902  if (enabled == false) {
903  this->a1 = a1;
904  this->a2 = a2;
905  enabled = true;
906  } else if (ref == "DEF") {
907  this->a1 = 0.0;
908  this->a2 = 0.0;
909  enabled = false;
910  } else {
911  // Otherwise, increment the values
912  this->a1 += a1;
913  this->a2 += a2;
914  enabled = true;
915  }
916  this->reference = ref;
917  }
918 
919  bool is_enabled() const {
920  return enabled;
921  };
922 
923  void to_json(rapidjson::Value& el, rapidjson::Document& doc) {
924  el.AddMember("type", "IdealHelmholtzEnthalpyEntropyOffset", doc.GetAllocator());
925  el.AddMember("a1", static_cast<double>(a1), doc.GetAllocator());
926  el.AddMember("a2", static_cast<double>(a2), doc.GetAllocator());
927  };
928  void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
929 };
930 
936 class IdealHelmholtzLogTau : public BaseHelmholtzTerm
937 {
938  private:
939  CoolPropDbl a1;
940  bool enabled;
941 
942  public:
944  IdealHelmholtzLogTau() : a1(_HUGE), enabled(false) {}
945 
946  // Constructor
947  IdealHelmholtzLogTau(CoolPropDbl a1) : a1(a1), enabled(true) {}
948 
949  bool is_enabled() const {
950  return enabled;
951  };
952 
953  void to_json(rapidjson::Value& el, rapidjson::Document& doc) {
954  el.AddMember("type", "IdealHelmholtzLogTau", doc.GetAllocator());
955  el.AddMember("a1", static_cast<double>(a1), doc.GetAllocator());
956  };
957  void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
958 };
959 
965 class IdealHelmholtzPower : public BaseHelmholtzTerm
966 {
967 
968  private:
969  std::vector<CoolPropDbl> n, t; // Use these variables internally
970  std::size_t N;
971  bool enabled;
972 
973  public:
974  IdealHelmholtzPower() : N(0), enabled(false){};
975  // Constructor
976  IdealHelmholtzPower(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& t) : n(n), t(t), N(n.size()), enabled(true){};
977 
978  bool is_enabled() const {
979  return enabled;
980  };
981 
982  void to_json(rapidjson::Value& el, rapidjson::Document& doc) {
983  el.AddMember("type", "IdealHelmholtzPower", doc.GetAllocator());
984  cpjson::set_long_double_array("n", n, el, doc);
985  cpjson::set_long_double_array("t", t, el, doc);
986  };
987  void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
988 };
989 
1036 class IdealHelmholtzPlanckEinsteinGeneralized : public BaseHelmholtzTerm
1037 {
1038 
1039  private:
1040  std::vector<CoolPropDbl> n, theta, c, d; // Use these variables internally
1041  std::size_t N;
1042  bool enabled;
1043 
1044  public:
1045  IdealHelmholtzPlanckEinsteinGeneralized() : N(0), enabled(false) {}
1046  // Constructor with std::vector instances
1047  IdealHelmholtzPlanckEinsteinGeneralized(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& theta,
1048  const std::vector<CoolPropDbl>& c, const std::vector<CoolPropDbl>& d)
1049  : n(n), theta(theta), c(c), d(d), N(n.size()), enabled(true) {}
1050 
1051  // Extend the vectors to allow for multiple instances feeding values to this function
1052  void extend(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& theta, const std::vector<CoolPropDbl>& c,
1053  const std::vector<CoolPropDbl>& d) {
1054  this->n.insert(this->n.end(), n.begin(), n.end());
1055  this->theta.insert(this->theta.end(), theta.begin(), theta.end());
1056  this->c.insert(this->c.end(), c.begin(), c.end());
1057  this->d.insert(this->d.end(), d.begin(), d.end());
1058  N += n.size();
1059  }
1060 
1061  bool is_enabled() const {
1062  return enabled;
1063  };
1064 
1065  void to_json(rapidjson::Value& el, rapidjson::Document& doc) {
1066  el.AddMember("type", "IdealHelmholtzPlanckEinsteinGeneralized", doc.GetAllocator());
1067  cpjson::set_long_double_array("n", n, el, doc);
1068  cpjson::set_long_double_array("theta", theta, el, doc);
1069  };
1070  void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
1071 };
1072 
1073 class IdealHelmholtzCP0Constant : public BaseHelmholtzTerm
1074 {
1075 
1076  private:
1077  double cp_over_R, Tc, T0, tau0; // Use these variables internally
1078  bool enabled;
1079 
1080  public:
1082  IdealHelmholtzCP0Constant() : cp_over_R(_HUGE), Tc(_HUGE), T0(_HUGE), tau0(_HUGE) {
1083  enabled = false;
1084  };
1085 
1087  IdealHelmholtzCP0Constant(CoolPropDbl cp_over_R, CoolPropDbl Tc, CoolPropDbl T0) : cp_over_R(cp_over_R), Tc(Tc), T0(T0) {
1088  enabled = true;
1089  tau0 = Tc / T0;
1090  };
1091 
1093  ~IdealHelmholtzCP0Constant(){};
1094 
1095  bool is_enabled() const {
1096  return enabled;
1097  };
1098 
1099  void to_json(rapidjson::Value& el, rapidjson::Document& doc) {
1100  el.AddMember("type", "IdealGasHelmholtzCP0Constant", doc.GetAllocator());
1101  el.AddMember("cp_over_R", cp_over_R, doc.GetAllocator());
1102  el.AddMember("Tc", Tc, doc.GetAllocator());
1103  el.AddMember("T0", T0, doc.GetAllocator());
1104  };
1105 
1106  void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
1107 };
1108 
1109 class IdealHelmholtzCP0PolyT : public BaseHelmholtzTerm
1110 {
1111  private:
1112  std::vector<CoolPropDbl> c, t;
1113  CoolPropDbl Tc, T0, tau0; // Use these variables internally
1114  std::size_t N;
1115  bool enabled;
1116 
1117  public:
1118  IdealHelmholtzCP0PolyT() : Tc(_HUGE), T0(_HUGE), tau0(_HUGE), N(0), enabled(false) {}
1119 
1121  IdealHelmholtzCP0PolyT(const std::vector<CoolPropDbl>& c, const std::vector<CoolPropDbl>& t, double Tc, double T0)
1122  : c(c), t(t), Tc(Tc), T0(T0), tau0(Tc / T0), N(c.size()), enabled(true) {
1123  assert(c.size() == t.size());
1124  }
1125 
1126  void extend(const std::vector<CoolPropDbl>& c, const std::vector<CoolPropDbl>& t) {
1127  this->c.insert(this->c.end(), c.begin(), c.end());
1128  this->t.insert(this->t.end(), t.begin(), t.end());
1129  N += c.size();
1130  }
1131 
1132  bool is_enabled() const {
1133  return enabled;
1134  };
1135 
1136  void to_json(rapidjson::Value& el, rapidjson::Document& doc);
1137  void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
1138 };
1142 class IdealHelmholtzGERG2004Sinh : public BaseHelmholtzTerm
1143 {
1144  private:
1145  std::vector<CoolPropDbl> n, theta;
1146  CoolPropDbl Tc, _Tr;
1147  std::size_t N;
1148  bool enabled;
1149 
1150  public:
1151  IdealHelmholtzGERG2004Sinh() : Tc(_HUGE), _Tr(_HUGE), N(0), enabled(false) {}
1152 
1154  IdealHelmholtzGERG2004Sinh(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& theta, double Tc)
1155  : n(n), theta(theta), Tc(Tc), _Tr(_HUGE), N(n.size()), enabled(true) {
1156  assert(n.size() == theta.size());
1157  }
1158 
1159  void extend(const std::vector<CoolPropDbl>& c, const std::vector<CoolPropDbl>& t) {
1160  this->n.insert(this->n.end(), n.begin(), n.end());
1161  this->theta.insert(this->theta.end(), theta.begin(), theta.end());
1162  N += c.size();
1163  }
1164  void set_Tred(CoolPropDbl Tr) {
1165  this->_Tr = Tr;
1166  }
1167 
1168  bool is_enabled() const {
1169  return enabled;
1170  };
1171  void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
1172 };
1173 
1174 class IdealHelmholtzGERG2004Cosh : public BaseHelmholtzTerm
1175 {
1176  private:
1177  std::vector<CoolPropDbl> n, theta;
1178  CoolPropDbl Tc, _Tr;
1179  std::size_t N;
1180  bool enabled;
1181 
1182  public:
1183  IdealHelmholtzGERG2004Cosh() : Tc(_HUGE), _Tr(_HUGE), N(0), enabled(false) {}
1184 
1186  IdealHelmholtzGERG2004Cosh(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& theta, double Tc)
1187  : n(n), theta(theta), Tc(Tc), _Tr(_HUGE), N(n.size()), enabled(true) {
1188  assert(n.size() == theta.size());
1189  }
1190 
1191  void extend(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& theta) {
1192  this->n.insert(this->n.end(), n.begin(), n.end());
1193  this->theta.insert(this->theta.end(), theta.begin(), theta.end());
1194  N += n.size();
1195  }
1196  void set_Tred(CoolPropDbl Tr) {
1197  this->_Tr = Tr;
1198  }
1199 
1200  bool is_enabled() const {
1201  return enabled;
1202  };
1203  void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
1204 };
1205 
1208 //\f[
1209 //\frac{c_p^0}{R_u} = A + B\left(\frac{C/T}{\sinh(C/T)}\right)^2 + D\left(\frac{E/T}{\cosh(E/T)}\right)^2
1210 //\f]
1211 //Second partial of ideal-gas Helmholtz energy given directly by specific heat (\f$\displaystyle\alpha_{\tau\tau}^0=-\frac{1}{\tau^2}\frac{c_p^0}{R_u} \f$) - this is obtained by real gas \f$c_p\f$ relationship, and killing off residual Helmholtz terms
1212 //\f[
1213 //\alpha^0_{\tau\tau} = -\frac{A}{\tau^2} - \frac{B}{\tau^2}\left(\frac{C/T}{\sinh(C/T)}\right)^2 - \frac{D}{\tau^2}\left(\frac{E/T}{\cosh(E/T)}\right)^2
1214 //\f]
1215 //or in terms of \f$ \tau \f$:
1216 //\f[
1217 //\alpha^0_{\tau\tau} = -\frac{A}{\tau^2} - \frac{BC^2}{T_c^2}\left(\frac{1}{\sinh(C\tau/T_c)}\right)^2 - \frac{DE^2}{T_c^2}\left(\frac{1}{\cosh(E\tau/T_c)}\right)^2
1218 //\f]
1219 //Third partial:
1220 //\f[
1221 //\alpha^0_{\tau\tau\tau} = 2\frac{A}{\tau^3} + 2\frac{BC^3}{T_c^3}\frac{\cosh(C\tau/T_c)}{\sinh^3(C\tau/T_c)} +2 \frac{DE^3}{T_c^3}\frac{\sinh(E\tau/T_c)}{\cosh^3(E\tau/T_c)}
1222 //\f]
1223 //Now coming back to the ideal gas Helmholtz energy definition:
1224 //\f[
1225 //\alpha^0 = -\tau\displaystyle\int_{\tau_0}^{\tau} \frac{1}{\tau^2}\frac{c_p^0}{R_u}d\tau+\displaystyle\int_{\tau_0}^{\tau} \frac{1}{\tau}\frac{c_p^0}{R_u}d\tau
1226 //\f]
1227 //Applying derivative
1228 //\f[
1229 //\alpha^0_{\tau} = -\displaystyle\int_{\tau_0}^{\tau} \frac{1}{\tau^2}\frac{c_p^0}{R_u}d\tau-\tau\frac{\partial}{\partial \tau}\left[\displaystyle\int_{\tau_0}^{\tau} \frac{1}{\tau^2}\frac{c_p^0}{R_u}d\tau \right]+\frac{\partial}{\partial \tau}\left[\displaystyle\int_{\tau_0}^{\tau} \frac{1}{\tau}\frac{c_p^0}{R_u}d\tau \right]
1230 //\f]
1231 //Fundamental theorem of calculus
1232 //\f[
1233 //\alpha^0_{\tau} = -\int_{\tau_0}^{\tau} \frac{1}{\tau^2}\frac{c_p^0}{R_u}d\tau-\tau \frac{1}{\tau^2}\frac{c_p^0}{R_u}d\tau+\frac{1}{\tau}\frac{c_p^0}{R_u}
1234 //\f]
1235 //Last two terms cancel, leaving
1236 //\f[
1237 //\alpha^0_{\tau} = -\int_{\tau_0}^{\tau} \frac{1}{\tau^2}\frac{c_p^0}{R_u}d\tau
1238 //\f]
1239 //Another derivative yields (from fundamental theorem of calculus)
1240 //\f[
1241 //\alpha^0_{\tau\tau} = - \frac{1}{\tau^2}\frac{c_p^0}{R_u}
1242 //\f]
1243 //
1244 //see also Jaeschke and Schley, 1995, (http://link.springer.com/article/10.1007%2FBF02083547#page-1)
1245 //*/
1247 //class IdealHelmholtzCP0AlyLee : public BaseHelmholtzTerm{
1248 //private:
1249 // std::vector<CoolPropDbl> c;
1250 // CoolPropDbl Tc, tau0, T0; // Use these variables internally
1251 // bool enabled;
1252 //public:
1253 // IdealHelmholtzCP0AlyLee(){enabled = false;};
1254 //
1255 // /// Constructor with std::vectors
1256 // IdealHelmholtzCP0AlyLee(const std::vector<CoolPropDbl> &c, double Tc, double T0)
1257 // :c(c), Tc(Tc), T0(T0)
1258 // {
1259 // tau0=Tc/T0;
1260 // enabled = true;
1261 // };
1262 //
1263 // /// Destructor
1264 // ~IdealHelmholtzCP0AlyLee(){};
1265 //
1266 // bool is_enabled() const {return enabled;};
1267 //
1268 // void to_json(rapidjson::Value &el, rapidjson::Document &doc);
1269 //
1270 //
1271 // /// The antiderivative given by \f$ \displaystyle\int \frac{1}{\tau^2}\frac{c_p^0}{R_u}d\tau \f$
1272 // /**
1273 // sympy code for this derivative:
1274 //
1275 // from sympy import *
1276 // a1,a2,a3,a4,a5,Tc,tau = symbols('a1,a2,a3,a4,a5,Tc,tau', real = True)
1277 // integrand = a1 + a2*(a3/Tc/sinh(a3*tau/Tc))**2 + a4*(a5/Tc/cosh(a5*tau/Tc))**2
1278 // integrand = integrand.rewrite(exp)
1279 // antideriv = trigsimp(integrate(integrand,tau))
1280 // display(antideriv)
1281 // print latex(antideriv)
1282 // print ccode(antideriv)
1283 //
1284 // \f[
1285 // \displaystyle\int \frac{1}{\tau^2}\frac{c_p^0}{R_u}d\tau = -\frac{a_0}{\tau}+\frac{2a_1a_2}{T_c\left[\exp\left(-\frac{2a_2\tau}{T_c}\right)-1\right]}+\frac{2a_3a_4}{T_c\left[\exp\left(-\frac{2a_4\tau}{T_c}\right)+1\right]}
1286 // \f]
1287 // */
1288 // CoolPropDbl anti_deriv_cp0_tau2(const CoolPropDbl &tau);
1289 //
1290 // /// The antiderivative given by \f$ \displaystyle\int \frac{1}{\tau}\frac{c_p^0}{R_u}d\tau \f$
1291 // /**
1292 // sympy code for this derivative:
1293 //
1294 // a_0,a_1,a_2,a_3,a_4,Tc,tau = symbols('a_0,a_1,a_2,a_3,a_4,Tc,tau', real = True)
1295 // integrand = a_0/tau + a_1/tau*(a_2*tau/Tc/sinh(a_2*tau/Tc))**2 + a_3/tau*(a_4*tau/Tc/cosh(a_4*tau/Tc))**2
1296 //
1297 // term2 = a_1/tau*(a_2*tau/Tc/sinh(a_2*tau/Tc))**2
1298 // term2 = term2.rewrite(exp) # Unpack the sinh to exp functions
1299 // antideriv2 = trigsimp(integrate(term2,tau))
1300 // display(antideriv2)
1301 // print latex(antideriv2)
1302 // print ccode(antideriv2)
1303 //
1304 // term3 = a_3/tau*(a_4*tau/Tc/cosh(a_4*tau/Tc))**2
1305 // term3 = term3.rewrite(exp) # Unpack the cosh to exp functions
1306 // antideriv3 = factor(trigsimp(integrate(term3,tau).rewrite(exp)))
1307 // display(antideriv3)
1308 // print latex(antideriv3)
1309 // print ccode(antideriv3)
1310 //
1311 // Can be broken into three parts (trick is to express \f$sinh\f$ and \f$cosh\f$ in terms of \f$exp\f$ function)
1312 //
1313 // Term 2:
1314 // \f[
1315 // \displaystyle\int \frac{a_1a_2^2}{T_c^2}\frac{\tau}{\sinh\left(\displaystyle\frac{a_2\tau}{T_c}\right)^2} d\tau = \frac{2 a_{1} a_{2} \tau}{- Tc + Tc e^{- \frac{2 a_{2}}{Tc} \tau}} + a_{1} \log{\left (-1 + e^{- \frac{2 a_{2}}{Tc} \tau} \right )} + \frac{2 a_{1}}{Tc} a_{2} \tau
1316 // \f]
1317 //
1318 // Term 3:
1319 // \f[
1320 // \displaystyle\int \frac{a_1a_2^2}{T_c^2}\frac{\tau}{\cosh\left(\displaystyle\frac{a_2\tau}{T_c}\right)^2} d\tau = - \frac{a_{3}}{Tc \left(e^{\frac{2 a_{4}}{Tc} \tau} + 1\right)} \left(Tc e^{\frac{2 a_{4}}{Tc} \tau} \log{\left (e^{\frac{2 a_{4}}{Tc} \tau} + 1 \right )} + Tc \log{\left (e^{\frac{2 a_{4}}{Tc} \tau} + 1 \right )} - 2 a_{4} \tau e^{\frac{2 a_{4}}{Tc} \tau}\right)
1321 // \f]
1322 // */
1323 // CoolPropDbl anti_deriv_cp0_tau(const CoolPropDbl &tau);
1324 //
1325 // CoolPropDbl base(const CoolPropDbl &tau, const CoolPropDbl &delta) throw();
1326 // CoolPropDbl dDelta(const CoolPropDbl &tau, const CoolPropDbl &delta) throw(){return 0.0;};
1327 // CoolPropDbl dTau(const CoolPropDbl &tau, const CoolPropDbl &delta) throw();
1328 // CoolPropDbl dDelta2(const CoolPropDbl &tau, const CoolPropDbl &delta) throw(){return 0.0;};
1329 // CoolPropDbl dDelta_dTau(const CoolPropDbl &tau, const CoolPropDbl &delta) throw(){return 0.0;};
1330 // CoolPropDbl dTau2(const CoolPropDbl &tau, const CoolPropDbl &delta) throw();
1331 // CoolPropDbl dDelta3(const CoolPropDbl &tau, const CoolPropDbl &delta) throw(){return 0.0;};
1332 // CoolPropDbl dDelta2_dTau(const CoolPropDbl &tau, const CoolPropDbl &delta) throw(){return 0.0;};
1333 // CoolPropDbl dDelta_dTau2(const CoolPropDbl &tau, const CoolPropDbl &delta) throw(){return 0.0;};
1334 // CoolPropDbl dTau3(const CoolPropDbl &tau, const CoolPropDbl &delta) throw();
1335 // CoolPropDbl dTau4(const CoolPropDbl &tau, const CoolPropDbl &delta) throw();
1336 //
1337 //};
1338 
1339 class IdealHelmholtzContainer : public BaseHelmholtzContainer
1340 {
1341  private:
1342  double _prefactor;
1343 
1344  public:
1345  IdealHelmholtzLead Lead;
1346  IdealHelmholtzEnthalpyEntropyOffset EnthalpyEntropyOffsetCore, EnthalpyEntropyOffset;
1347  IdealHelmholtzLogTau LogTau;
1348  IdealHelmholtzPower Power;
1349  IdealHelmholtzPlanckEinsteinGeneralized PlanckEinstein;
1350 
1351  IdealHelmholtzCP0Constant CP0Constant;
1352  IdealHelmholtzCP0PolyT CP0PolyT;
1353  IdealHelmholtzGERG2004Cosh GERG2004Cosh;
1354  IdealHelmholtzGERG2004Sinh GERG2004Sinh;
1355 
1356  IdealHelmholtzContainer() : _prefactor(1.0){};
1357 
1358  void set_prefactor(double prefactor) {
1359  _prefactor = prefactor;
1360  }
1361 
1362  void set_Tred(double T_red) {
1363  GERG2004Cosh.set_Tred(T_red);
1364  GERG2004Sinh.set_Tred(T_red);
1365  }
1366 
1367  void empty_the_EOS() {
1368  Lead = IdealHelmholtzLead();
1369  EnthalpyEntropyOffsetCore = IdealHelmholtzEnthalpyEntropyOffset();
1370  EnthalpyEntropyOffset = IdealHelmholtzEnthalpyEntropyOffset();
1371  LogTau = IdealHelmholtzLogTau();
1372  Power = IdealHelmholtzPower();
1373  PlanckEinstein = IdealHelmholtzPlanckEinsteinGeneralized();
1374  CP0Constant = IdealHelmholtzCP0Constant();
1375  CP0PolyT = IdealHelmholtzCP0PolyT();
1376  GERG2004Cosh = IdealHelmholtzGERG2004Cosh();
1377  GERG2004Sinh = IdealHelmholtzGERG2004Sinh();
1378  };
1379 
1380  HelmholtzDerivatives all(const CoolPropDbl tau, const CoolPropDbl delta, bool cache_values = false) {
1381  HelmholtzDerivatives derivs; // zeros out the elements
1382  Lead.all(tau, delta, derivs);
1383  EnthalpyEntropyOffsetCore.all(tau, delta, derivs);
1384  EnthalpyEntropyOffset.all(tau, delta, derivs);
1385  LogTau.all(tau, delta, derivs);
1386  Power.all(tau, delta, derivs);
1387  PlanckEinstein.all(tau, delta, derivs);
1388  CP0Constant.all(tau, delta, derivs);
1389  CP0PolyT.all(tau, delta, derivs);
1390  GERG2004Cosh.all(tau, delta, derivs);
1391  GERG2004Sinh.all(tau, delta, derivs);
1392 
1393  if (cache_values) {
1394  _base = derivs.alphar * _prefactor;
1395  _dDelta = derivs.dalphar_ddelta * _prefactor;
1396  _dTau = derivs.dalphar_dtau * _prefactor;
1397  _dDelta2 = derivs.d2alphar_ddelta2 * _prefactor;
1398  _dTau2 = derivs.d2alphar_dtau2 * _prefactor;
1399  _dDelta_dTau = derivs.d2alphar_ddelta_dtau * _prefactor;
1400  _dDelta3 = derivs.d3alphar_ddelta3 * _prefactor;
1401  _dTau3 = derivs.d3alphar_dtau3 * _prefactor;
1402  _dDelta2_dTau = derivs.d3alphar_ddelta2_dtau * _prefactor;
1403  _dDelta_dTau2 = derivs.d3alphar_ddelta_dtau2 * _prefactor;
1404  }
1405  return derivs * _prefactor;
1406  };
1407 };
1408 }; /* namespace CoolProp */
1409 
1410 #endif
Generated on Wed Nov 29 2023 14:26:23 for CoolProp by doxygen 1.9.1