Loading [MathJax]/extensions/TeX/AMSmath.js
CoolProp  6.7.1dev
An open-source fluid property and humid air property database
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 #if ENABLE_CATCH
14 # include "MultiComplex/MultiComplex.hpp"
15 #endif
16 
17 namespace CoolProp {
18 
19 // #############################################################################
20 // #############################################################################
21 // #############################################################################
22 // RESIDUAL TERMS
23 // #############################################################################
24 // #############################################################################
25 // #############################################################################
26 
27 #define LIST_OF_DERIVATIVE_VARIABLES \
28  X(alphar) \
29  X(dalphar_ddelta) \
30  X(dalphar_dtau) \
31  X(d2alphar_ddelta2) \
32  X(d2alphar_dtau2) \
33  X(d2alphar_ddelta_dtau) \
34  X(d3alphar_ddelta3) \
35  X(d3alphar_ddelta_dtau2) \
36  X(d3alphar_ddelta2_dtau) \
37  X(d3alphar_dtau3) \
38  X(d4alphar_ddelta4) \
39  X(d4alphar_ddelta3_dtau) \
40  X(d4alphar_ddelta2_dtau2) \
41  X(d4alphar_ddelta_dtau3) \
42  X(d4alphar_dtau4) \
43  X(delta_x_dalphar_ddelta) \
44  X(tau_x_dalphar_dtau) \
45  X(delta2_x_d2alphar_ddelta2) \
46  X(deltatau_x_d2alphar_ddelta_dtau) \
47  X(tau2_x_d2alphar_dtau2)
48 
49 struct HelmholtzDerivatives
50 {
51 #define X(name) CoolPropDbl name;
52  LIST_OF_DERIVATIVE_VARIABLES
53 #undef X
54  CoolPropDbl tau, delta, T_red, rhomolar_red;
55 
56  void reset(CoolPropDbl v) {
57 #define X(name) name = v;
58  LIST_OF_DERIVATIVE_VARIABLES
59 #undef X
60  }
61  HelmholtzDerivatives operator+(const HelmholtzDerivatives& other) const {
62  HelmholtzDerivatives _new;
63 #define X(name) _new.name = name + other.name;
64  LIST_OF_DERIVATIVE_VARIABLES
65 #undef X
66  return _new;
67  }
68  HelmholtzDerivatives operator*(const CoolPropDbl& other) const {
69  HelmholtzDerivatives _new;
70 #define X(name) _new.name = name * other;
71  LIST_OF_DERIVATIVE_VARIABLES
72 #undef X
73  return _new;
74  }
75  HelmholtzDerivatives() {
76  reset(0.0);
77  T_red = _HUGE;
78  rhomolar_red = _HUGE;
79  };
81  double get(std::size_t itau, std::size_t idelta) {
82  if (itau == 0) {
83  if (idelta == 0) {
84  return alphar;
85  } else if (idelta == 1) {
86  return dalphar_ddelta;
87  } else if (idelta == 2) {
88  return d2alphar_ddelta2;
89  } else if (idelta == 3) {
90  return d3alphar_ddelta3;
91  } else if (idelta == 4) {
92  return d4alphar_ddelta4;
93  } else {
94  throw ValueError();
95  }
96  } else if (itau == 1) {
97  if (idelta == 0) {
98  return dalphar_dtau;
99  } else if (idelta == 1) {
100  return d2alphar_ddelta_dtau;
101  } else if (idelta == 2) {
102  return d3alphar_ddelta2_dtau;
103  } else if (idelta == 3) {
104  return d4alphar_ddelta3_dtau;
105  } else {
106  throw ValueError();
107  }
108  } else if (itau == 2) {
109  if (idelta == 0) {
110  return d2alphar_dtau2;
111  } else if (idelta == 1) {
112  return d3alphar_ddelta_dtau2;
113  } else if (idelta == 2) {
114  return d4alphar_ddelta2_dtau2;
115  } else {
116  throw ValueError();
117  }
118  } else if (itau == 3) {
119  if (idelta == 0) {
120  return d3alphar_dtau3;
121  } else if (idelta == 1) {
122  return d4alphar_ddelta_dtau3;
123  } else {
124  throw ValueError();
125  }
126  } else if (itau == 4) {
127  if (idelta == 0) {
128  return d4alphar_dtau4;
129  } else {
130  throw ValueError();
131  }
132  } else {
133  throw ValueError();
134  }
135  }
136 };
137 #undef LIST_OF_DERIVATIVE_VARIABLES
138 
140 
165 class BaseHelmholtzTerm
166 {
167  public:
168  BaseHelmholtzTerm(){};
169  virtual ~BaseHelmholtzTerm(){};
170 
172 
175  virtual CoolPropDbl base(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
176  HelmholtzDerivatives deriv;
177  all(tau, delta, deriv);
178  return deriv.alphar;
179  };
181 
184  virtual CoolPropDbl dTau(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
185  HelmholtzDerivatives deriv;
186  all(tau, delta, deriv);
187  return deriv.dalphar_dtau;
188  };
190 
193  virtual CoolPropDbl dTau2(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
194  HelmholtzDerivatives deriv;
195  all(tau, delta, deriv);
196  return deriv.d2alphar_dtau2;
197  };
199 
202  virtual CoolPropDbl dDelta_dTau(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
203  HelmholtzDerivatives deriv;
204  all(tau, delta, deriv);
205  return deriv.d2alphar_ddelta_dtau;
206  };
208 
211  virtual CoolPropDbl dDelta(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
212  HelmholtzDerivatives deriv;
213  all(tau, delta, deriv);
214  return deriv.dalphar_ddelta;
215  };
217 
220  virtual CoolPropDbl dDelta2(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
221  HelmholtzDerivatives deriv;
222  all(tau, delta, deriv);
223  return deriv.d2alphar_ddelta2;
224  };
226 
229  virtual CoolPropDbl dDelta2_dTau(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
230  HelmholtzDerivatives deriv;
231  all(tau, delta, deriv);
232  return deriv.d3alphar_ddelta2_dtau;
233  };
235 
238  virtual CoolPropDbl dDelta_dTau2(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
239  HelmholtzDerivatives deriv;
240  all(tau, delta, deriv);
241  return deriv.d3alphar_ddelta_dtau2;
242  };
244 
247  virtual CoolPropDbl dTau3(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
248  HelmholtzDerivatives deriv;
249  all(tau, delta, deriv);
250  return deriv.d3alphar_dtau3;
251  };
253 
256  virtual CoolPropDbl dDelta3(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
257  HelmholtzDerivatives deriv;
258  all(tau, delta, deriv);
259  return deriv.d3alphar_ddelta3;
260  };
262 
265  virtual CoolPropDbl dTau4(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
266  HelmholtzDerivatives deriv;
267  all(tau, delta, deriv);
268  return deriv.d4alphar_dtau4;
269  };
270  virtual CoolPropDbl dDelta_dTau3(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
271  HelmholtzDerivatives deriv;
272  all(tau, delta, deriv);
273  return deriv.d4alphar_ddelta_dtau3;
274  };
275  virtual CoolPropDbl dDelta2_dTau2(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
276  HelmholtzDerivatives deriv;
277  all(tau, delta, deriv);
278  return deriv.d4alphar_ddelta2_dtau2;
279  };
280  virtual CoolPropDbl dDelta3_dTau(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
281  HelmholtzDerivatives deriv;
282  all(tau, delta, deriv);
283  return deriv.d4alphar_ddelta3_dtau;
284  };
285  virtual CoolPropDbl dDelta4(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
286  HelmholtzDerivatives deriv;
287  all(tau, delta, deriv);
288  return deriv.d4alphar_ddelta4;
289  };
290 
291  virtual void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw() = 0;
292  #if ENABLE_CATCH
293  virtual mcx::MultiComplex<double> one_mcx(const mcx::MultiComplex<double>& tau, const mcx::MultiComplex<double>& delta) const {
294  throw CoolProp::NotImplementedError("The mcx derivative function was not implemented");
295  }
296  #endif
297 };
298 
299 struct ResidualHelmholtzGeneralizedExponentialElement
300 {
302  CoolPropDbl n, d, t;
305  CoolPropDbl c, l_double, omega, m_double, eta1, epsilon1, eta2, epsilon2, beta1, gamma1, beta2, gamma2;
307  int l_int, m_int;
309  bool l_is_int, m_is_int;
310 
311  ResidualHelmholtzGeneralizedExponentialElement() {
312  n = 0;
313  d = 0;
314  t = 0;
315  c = 0;
316  l_double = 0;
317  omega = 0;
318  m_double = 0;
319  eta1 = 0;
320  epsilon1 = 0;
321  eta2 = 0;
322  epsilon2 = 0;
323  beta1 = 0;
324  gamma1 = 0;
325  beta2 = 0;
326  gamma2 = 0;
327  l_int = 0;
328  m_int = 0;
329  l_is_int = false;
330  m_is_int = true;
331  };
332 };
341 class ResidualHelmholtzGeneralizedExponential : public BaseHelmholtzTerm
342 {
343 
344  public:
345  bool delta_li_in_u, tau_mi_in_u, eta1_in_u, eta2_in_u, beta1_in_u, beta2_in_u, finished;
346  std::vector<CoolPropDbl> s;
347  std::size_t N;
348 
349  // These variables are for the exp(u) part
350  // 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
351  std::vector<double> n, d, t, c, l_double, omega, m_double, eta1, epsilon1, eta2, epsilon2, beta1, gamma1, beta2, gamma2;
352  // 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)
353  std::vector<int> l_int, m_int;
354 
355  //Eigen::ArrayXd uE, du_ddeltaE, du_dtauE, d2u_ddelta2E, d2u_dtau2E, d3u_ddelta3E, d3u_dtau3E;
356 
357  std::vector<ResidualHelmholtzGeneralizedExponentialElement> elements;
358  // Default Constructor
359  ResidualHelmholtzGeneralizedExponential()
360  : 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){};
366  void add_Power(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& d, const std::vector<CoolPropDbl>& t,
367  const std::vector<CoolPropDbl>& l) {
368  for (std::size_t i = 0; i < n.size(); ++i) {
369  ResidualHelmholtzGeneralizedExponentialElement el;
370  el.n = n[i];
371  el.d = d[i];
372  el.t = t[i];
373  el.l_double = l[i];
374  el.l_int = (int)el.l_double;
375  if (el.l_double > 0)
376  el.c = 1.0;
377  else
378  el.c = 0.0;
379  elements.push_back(el);
380  }
381  delta_li_in_u = true;
382  };
388  void add_Exponential(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& d, const std::vector<CoolPropDbl>& t,
389  const std::vector<CoolPropDbl>& g, const std::vector<CoolPropDbl>& l) {
390  for (std::size_t i = 0; i < n.size(); ++i) {
391  ResidualHelmholtzGeneralizedExponentialElement el;
392  el.n = n[i];
393  el.d = d[i];
394  el.t = t[i];
395  el.c = g[i];
396  el.l_double = l[i];
397  el.l_int = (int)el.l_double;
398  elements.push_back(el);
399  }
400  delta_li_in_u = true;
401  }
407  void add_Gaussian(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& d, const std::vector<CoolPropDbl>& t,
408  const std::vector<CoolPropDbl>& eta, const std::vector<CoolPropDbl>& epsilon, const std::vector<CoolPropDbl>& beta,
409  const std::vector<CoolPropDbl>& gamma) {
410  for (std::size_t i = 0; i < n.size(); ++i) {
411  ResidualHelmholtzGeneralizedExponentialElement el;
412  el.n = n[i];
413  el.d = d[i];
414  el.t = t[i];
415  el.eta2 = eta[i];
416  el.epsilon2 = epsilon[i];
417  el.beta2 = beta[i];
418  el.gamma2 = gamma[i];
419  elements.push_back(el);
420  }
421  eta2_in_u = true;
422  beta2_in_u = true;
423  };
429  void add_GERG2008Gaussian(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& d, const std::vector<CoolPropDbl>& t,
430  const std::vector<CoolPropDbl>& eta, const std::vector<CoolPropDbl>& epsilon, const std::vector<CoolPropDbl>& beta,
431  const std::vector<CoolPropDbl>& gamma) {
432  for (std::size_t i = 0; i < n.size(); ++i) {
433  ResidualHelmholtzGeneralizedExponentialElement el;
434  el.n = n[i];
435  el.d = d[i];
436  el.t = t[i];
437  el.eta2 = eta[i];
438  el.epsilon2 = epsilon[i];
439  el.eta1 = beta[i];
440  el.epsilon1 = gamma[i];
441  elements.push_back(el);
442  }
443  eta2_in_u = true;
444  eta1_in_u = true;
445  };
451  void add_Lemmon2005(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& d, const std::vector<CoolPropDbl>& t,
452  const std::vector<CoolPropDbl>& l, const std::vector<CoolPropDbl>& m) {
453  for (std::size_t i = 0; i < n.size(); ++i) {
454  ResidualHelmholtzGeneralizedExponentialElement el;
455  el.n = n[i];
456  el.d = d[i];
457  el.t = t[i];
458  el.c = 1.0;
459  el.omega = 1.0;
460  el.l_double = l[i];
461  el.m_double = m[i];
462  el.l_int = (int)el.l_double;
463  el.m_int = (int)el.m_double;
464  elements.push_back(el);
465  }
466  delta_li_in_u = true;
467  tau_mi_in_u = true;
468  };
469 
470  void finish() {
471  n.resize(elements.size());
472  d.resize(elements.size());
473  t.resize(elements.size());
474  c.resize(elements.size());
475  omega.resize(elements.size());
476  l_double.resize(elements.size());
477  l_int.resize(elements.size());
478  m_double.resize(elements.size());
479  m_int.resize(elements.size());
480  epsilon2.resize(elements.size());
481  eta2.resize(elements.size());
482  gamma2.resize(elements.size());
483  beta2.resize(elements.size());
484 
485  for (std::size_t i = 0; i < elements.size(); ++i) {
486  n[i] = elements[i].n;
487  d[i] = elements[i].d;
488  t[i] = elements[i].t;
489  c[i] = elements[i].c;
490  omega[i] = elements[i].omega;
491  l_double[i] = elements[i].l_double;
492  l_int[i] = elements[i].l_int;
493  m_double[i] = elements[i].m_double;
494  m_int[i] = elements[i].m_int;
495  epsilon2[i] = elements[i].epsilon2;
496  eta2[i] = elements[i].eta2;
497  gamma2[i] = elements[i].gamma2;
498  beta2[i] = elements[i].beta2;
499 
500  // See if l is an integer, and store a flag if it is
501  elements[i].l_is_int = (std::abs(static_cast<long>(elements[i].l_double) - elements[i].l_double) < 1e-14);
502  }
503  // uE.resize(elements.size());
504  // du_ddeltaE.resize(elements.size());
505  // du_dtauE.resize(elements.size());
506  // d2u_ddelta2E.resize(elements.size());
507  // d2u_dtau2E.resize(elements.size());
508  // d3u_ddelta3E.resize(elements.size());
509  // d3u_dtau3E.resize(elements.size());
510 
511  finished = true;
512  };
513 
514  void to_json(rapidjson::Value& el, rapidjson::Document& doc);
515 
516  void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
517  //void allEigen(const CoolPropDbl &tau, const CoolPropDbl &delta, HelmholtzDerivatives &derivs) throw();
518 
519  #if ENABLE_CATCH
520  virtual mcx::MultiComplex<double> one_mcx(const mcx::MultiComplex<double>& tau, const mcx::MultiComplex<double>& delta) const override;
521  #endif
522 };
523 
524 struct ResidualHelmholtzNonAnalyticElement
525 {
526  CoolPropDbl n, a, b, beta, A, B, C, D;
527 };
528 class ResidualHelmholtzNonAnalytic : public BaseHelmholtzTerm
529 {
530 
531  public:
532  std::size_t N;
533  std::vector<CoolPropDbl> s;
534  std::vector<ResidualHelmholtzNonAnalyticElement> elements;
536  ResidualHelmholtzNonAnalytic() {
537  N = 0;
538  };
540  ~ResidualHelmholtzNonAnalytic(){};
542  ResidualHelmholtzNonAnalytic(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& a, const std::vector<CoolPropDbl>& b,
543  const std::vector<CoolPropDbl>& beta, const std::vector<CoolPropDbl>& A, const std::vector<CoolPropDbl>& B,
544  const std::vector<CoolPropDbl>& C, const std::vector<CoolPropDbl>& D) {
545  N = n.size();
546  s.resize(N);
547  for (std::size_t i = 0; i < n.size(); ++i) {
548  ResidualHelmholtzNonAnalyticElement el;
549  el.n = n[i];
550  el.a = a[i];
551  el.b = b[i];
552  el.beta = beta[i];
553  el.A = A[i];
554  el.B = B[i];
555  el.C = C[i];
556  el.D = D[i];
557  elements.push_back(el);
558  }
559  };
560  void to_json(rapidjson::Value& el, rapidjson::Document& doc);
561  void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
562 #if ENABLE_CATCH
563  virtual mcx::MultiComplex<double> one_mcx(const mcx::MultiComplex<double>& tau, const mcx::MultiComplex<double>& delta) const override;
564 #endif
565 };
566 
567 class ResidualHelmholtzGeneralizedCubic : public BaseHelmholtzTerm
568 {
569  protected:
570  shared_ptr<AbstractCubic> m_abstractcubic;
571  std::vector<double> z;
572  public:
573  bool enabled;
574 
576  ResidualHelmholtzGeneralizedCubic() {
577  enabled = false;
578  };
580  ResidualHelmholtzGeneralizedCubic(shared_ptr<AbstractCubic>& ac) : m_abstractcubic(ac) {
581  enabled = true;
582  z = std::vector<double>(1, 1); // Init the vector to [1.0]
583  };
584 
585  void to_json(rapidjson::Value& el, rapidjson::Document& doc);
586  void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
587 };
588 
589 class ResidualHelmholtzGaoB : public BaseHelmholtzTerm
590 {
591  protected:
592  std::vector<double> n, t, d, eta, beta, gamma, epsilon, b;
593 
594  public:
595  bool enabled;
596 
598  ResidualHelmholtzGaoB() {
599  enabled = false;
600  };
601 
603  ResidualHelmholtzGaoB(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& t, const std::vector<CoolPropDbl>& d,
604  const std::vector<CoolPropDbl>& eta, const std::vector<CoolPropDbl>& beta, const std::vector<CoolPropDbl>& gamma,
605  const std::vector<CoolPropDbl>& epsilon, const std::vector<CoolPropDbl>& b)
606  : n(n), t(t), d(d), eta(eta), beta(beta), gamma(gamma), epsilon(epsilon), b(b) {
607  enabled = true;
608  };
609 
610  void to_json(rapidjson::Value& el, rapidjson::Document& doc);
611  void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
612 
613  #if ENABLE_CATCH
614  virtual mcx::MultiComplex<double> one_mcx(const mcx::MultiComplex<double>& tau, const mcx::MultiComplex<double>& delta) const override;
615  #endif
616 };
617 
619 class ResidualHelmholtzXiangDeiters : public BaseHelmholtzTerm
620 {
621 
622  public:
623  bool enabled;
624  ResidualHelmholtzGeneralizedExponential phi0, phi1, phi2;
625  CoolPropDbl Tc, pc, rhomolarc, acentric, R, theta;
627  ResidualHelmholtzXiangDeiters() : Tc(_HUGE), pc(_HUGE), rhomolarc(_HUGE), acentric(_HUGE), R(_HUGE), theta(_HUGE) {
628  enabled = false;
629  };
631  ResidualHelmholtzXiangDeiters(const CoolPropDbl Tc, const CoolPropDbl pc, const CoolPropDbl rhomolarc, const CoolPropDbl acentric,
632  const CoolPropDbl R);
633  void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
634 #if ENABLE_CATCH
635  virtual mcx::MultiComplex<double> one_mcx(const mcx::MultiComplex<double>& tau, const mcx::MultiComplex<double>& delta) const override;
636 #endif
637 };
638 
639 class ResidualHelmholtzSAFTAssociating : public BaseHelmholtzTerm
640 {
641 
642  protected:
643  double a, m, epsilonbar, vbarn, kappabar;
644 
645  CoolPropDbl Deltabar(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
646  CoolPropDbl dDeltabar_ddelta__consttau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
647  CoolPropDbl d2Deltabar_ddelta2__consttau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
648  CoolPropDbl dDeltabar_dtau__constdelta(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
649  CoolPropDbl d2Deltabar_dtau2__constdelta(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
650  CoolPropDbl d2Deltabar_ddelta_dtau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
651  CoolPropDbl d3Deltabar_dtau3__constdelta(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
652  CoolPropDbl d3Deltabar_ddelta_dtau2(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
653  CoolPropDbl d3Deltabar_ddelta3__consttau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
654  CoolPropDbl d3Deltabar_ddelta2_dtau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
655 
656  CoolPropDbl X(const CoolPropDbl& delta, const CoolPropDbl& Deltabar) const;
657  CoolPropDbl dX_dDeltabar__constdelta(const CoolPropDbl& delta, const CoolPropDbl& Deltabar) const;
658  CoolPropDbl dX_ddelta__constDeltabar(const CoolPropDbl& delta, const CoolPropDbl& Deltabar) const;
659  CoolPropDbl dX_dtau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
660  CoolPropDbl dX_ddelta(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
661  CoolPropDbl d2X_dtau2(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
662  CoolPropDbl d2X_ddeltadtau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
663  CoolPropDbl d2X_ddelta2(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
664 
665  CoolPropDbl d3X_dtau3(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
666  CoolPropDbl d3X_ddelta3(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
667  CoolPropDbl d3X_ddeltadtau2(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
668  CoolPropDbl d3X_ddelta2dtau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
669 
670  CoolPropDbl g(const CoolPropDbl& eta) const;
671  CoolPropDbl dg_deta(const CoolPropDbl& eta) const;
672  CoolPropDbl d2g_deta2(const CoolPropDbl& eta) const;
673  CoolPropDbl d3g_deta3(const CoolPropDbl& eta) const;
674  CoolPropDbl eta(const CoolPropDbl& delta) const;
675 
676  public:
678  ResidualHelmholtzSAFTAssociating() : a(_HUGE), m(_HUGE), epsilonbar(_HUGE), vbarn(_HUGE), kappabar(_HUGE) {
679  disabled = true;
680  };
681 
682  // Constructor
683  ResidualHelmholtzSAFTAssociating(double a, double m, double epsilonbar, double vbarn, double kappabar)
684  : a(a), m(m), epsilonbar(epsilonbar), vbarn(vbarn), kappabar(kappabar) {
685  disabled = false;
686  };
687 
688  bool disabled;
689 
690  //Destructor. No Implementation
691  ~ResidualHelmholtzSAFTAssociating(){};
692 
693  void to_json(rapidjson::Value& el, rapidjson::Document& doc);
694 
695  CoolPropDbl dTau4(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
696  return 1e99;
697  };
698  CoolPropDbl dDelta_dTau3(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
699  return 1e99;
700  };
701  CoolPropDbl dDelta2_dTau2(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
702  return 1e99;
703  };
704  CoolPropDbl dDelta3_dTau(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
705  return 1e99;
706  };
707  CoolPropDbl dDelta4(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
708  return 1e99;
709  };
710 
711  void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& deriv) throw();
712 };
713 
714 class BaseHelmholtzContainer
715 {
716  protected:
717  CachedElement _base, _dDelta, _dTau, _dDelta2, _dTau2, _dDelta_dTau, _dDelta3, _dDelta2_dTau, _dDelta_dTau2, _dTau3;
718  CachedElement _dDelta4, _dDelta3_dTau, _dDelta2_dTau2, _dDelta_dTau3, _dTau4;
719 
720  public:
721  void clear() {
722  _base.clear();
723  _dDelta.clear();
724  _dTau.clear();
725  _dDelta2.clear();
726  _dTau2.clear();
727  _dDelta_dTau.clear();
728  _dDelta3.clear();
729  _dTau3.clear();
730  _dDelta2_dTau.clear();
731  _dDelta_dTau2.clear();
732  _dDelta4.clear();
733  _dDelta3_dTau.clear();
734  _dDelta2_dTau2.clear();
735  _dDelta_dTau3.clear();
736  _dTau4.clear();
737  };
738 
739  virtual void empty_the_EOS() = 0;
740  virtual HelmholtzDerivatives all(const CoolPropDbl tau, const CoolPropDbl delta, bool cache_values) = 0;
741 
742  CoolPropDbl base(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
743  if (!_base || dont_use_cache)
744  return all(tau, delta, false).alphar;
745  else
746  return _base;
747  };
748  CoolPropDbl dDelta(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
749  if (!_dDelta || dont_use_cache)
750  return all(tau, delta, false).dalphar_ddelta;
751  else
752  return _dDelta;
753  };
754  CoolPropDbl dTau(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
755  if (!_dTau || dont_use_cache)
756  return all(tau, delta, false).dalphar_dtau;
757  else
758  return _dTau;
759  };
760  CoolPropDbl dDelta2(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
761  if (!_dDelta2 || dont_use_cache)
762  return all(tau, delta, false).d2alphar_ddelta2;
763  else
764  return _dDelta2;
765  };
766  CoolPropDbl dDelta_dTau(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
767  if (!_dDelta_dTau || dont_use_cache)
768  return all(tau, delta, false).d2alphar_ddelta_dtau;
769  else
770  return _dDelta_dTau;
771  };
772  CoolPropDbl dTau2(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
773  if (!_dTau2 || dont_use_cache)
774  return all(tau, delta, false).d2alphar_dtau2;
775  else
776  return _dTau2;
777  };
778  CoolPropDbl dDelta3(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
779  if (!_dDelta3 || dont_use_cache)
780  return all(tau, delta, false).d3alphar_ddelta3;
781  else
782  return _dDelta3;
783  };
784  CoolPropDbl dDelta2_dTau(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
785  if (!_dDelta2_dTau || dont_use_cache)
786  return all(tau, delta, false).d3alphar_ddelta2_dtau;
787  else
788  return _dDelta2_dTau;
789  };
790  CoolPropDbl dDelta_dTau2(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
791  if (!_dDelta_dTau2 || dont_use_cache)
792  return all(tau, delta, false).d3alphar_ddelta_dtau2;
793  else
794  return _dDelta_dTau2;
795  };
796  CoolPropDbl dTau3(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
797  if (!_dTau3 || dont_use_cache)
798  return all(tau, delta, false).d3alphar_dtau3;
799  else
800  return _dTau3;
801  };
802  CoolPropDbl dDelta4(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
803  return all(tau, delta, false).d4alphar_ddelta4;
804  };
805  CoolPropDbl dDelta3_dTau(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
806  return all(tau, delta, false).d4alphar_ddelta3_dtau;
807  };
808  CoolPropDbl dDelta2_dTau2(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
809  return all(tau, delta, false).d4alphar_ddelta2_dtau2;
810  };
811  CoolPropDbl dDelta_dTau3(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
812  return all(tau, delta, false).d4alphar_ddelta_dtau3;
813  };
814  CoolPropDbl dTau4(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
815  return all(tau, delta, false).d4alphar_dtau4;
816  };
817 };
818 
819 class ResidualHelmholtzContainer : public BaseHelmholtzContainer
820 {
821  public:
822  ResidualHelmholtzNonAnalytic NonAnalytic;
823  ResidualHelmholtzSAFTAssociating SAFT;
824  ResidualHelmholtzGeneralizedExponential GenExp;
825  ResidualHelmholtzGeneralizedCubic cubic;
826  ResidualHelmholtzXiangDeiters XiangDeiters;
827  ResidualHelmholtzGaoB GaoB;
828 
829  void empty_the_EOS() {
830  NonAnalytic = ResidualHelmholtzNonAnalytic();
831  SAFT = ResidualHelmholtzSAFTAssociating();
832  GenExp = ResidualHelmholtzGeneralizedExponential();
833  cubic = ResidualHelmholtzGeneralizedCubic();
834  XiangDeiters = ResidualHelmholtzXiangDeiters();
835  GaoB = ResidualHelmholtzGaoB();
836  };
837 
838  HelmholtzDerivatives all(const CoolPropDbl tau, const CoolPropDbl delta, bool cache_values = false) {
839  HelmholtzDerivatives derivs; // zeros out the elements
840  GenExp.all(tau, delta, derivs);
841  NonAnalytic.all(tau, delta, derivs);
842  SAFT.all(tau, delta, derivs);
843  cubic.all(tau, delta, derivs);
844  XiangDeiters.all(tau, delta, derivs);
845  GaoB.all(tau, delta, derivs);
846  if (cache_values) {
847  _base = derivs.alphar;
848  _dDelta = derivs.dalphar_ddelta;
849  _dTau = derivs.dalphar_dtau;
850  _dDelta2 = derivs.d2alphar_ddelta2;
851  _dTau2 = derivs.d2alphar_dtau2;
852  _dDelta_dTau = derivs.d2alphar_ddelta_dtau;
853  _dDelta3 = derivs.d3alphar_ddelta3;
854  _dTau3 = derivs.d3alphar_dtau3;
855  _dDelta2_dTau = derivs.d3alphar_ddelta2_dtau;
856  _dDelta_dTau2 = derivs.d3alphar_ddelta_dtau2;
857  }
858  return derivs;
859  };
860 };
861 
862 // #############################################################################
863 // #############################################################################
864 // #############################################################################
865 // IDEAL GAS TERMS
866 // #############################################################################
867 // #############################################################################
868 // #############################################################################
869 
871 
876 class IdealHelmholtzLead : public BaseHelmholtzTerm
877 {
878 
879  private:
880  CoolPropDbl a1, a2;
881  bool enabled;
882 
883  public:
884  // Default constructor
885  IdealHelmholtzLead() : a1(_HUGE), a2(_HUGE), enabled(false) {}
886 
887  // Constructor
888  IdealHelmholtzLead(CoolPropDbl a1, CoolPropDbl a2) : a1(a1), a2(a2), enabled(true) {}
889 
890  bool is_enabled() const {
891  return enabled;
892  }
893 
894  void to_json(rapidjson::Value& el, rapidjson::Document& doc) {
895  el.AddMember("type", "IdealHelmholtzLead", doc.GetAllocator());
896  el.AddMember("a1", static_cast<double>(a1), doc.GetAllocator());
897  el.AddMember("a2", static_cast<double>(a2), doc.GetAllocator());
898  };
899 
900  void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
901 };
902 
904 
909 class IdealHelmholtzEnthalpyEntropyOffset : public BaseHelmholtzTerm
910 {
911  private:
912  CoolPropDbl a1, a2; // Use these variables internally
913  std::string reference;
914  bool enabled;
915 
916  public:
917  IdealHelmholtzEnthalpyEntropyOffset() : a1(_HUGE), a2(_HUGE), enabled(false) {}
918 
919  // Constructor
920  IdealHelmholtzEnthalpyEntropyOffset(CoolPropDbl a1, CoolPropDbl a2, const std::string& ref) : a1(a1), a2(a2), reference(ref), enabled(true) {}
921 
922  // Set the values in the class
923  void set(CoolPropDbl a1, CoolPropDbl a2, const std::string& ref) {
924  // If it doesn't already exist, just set the values
925  if (enabled == false) {
926  this->a1 = a1;
927  this->a2 = a2;
928  enabled = true;
929  } else if (ref == "DEF") {
930  this->a1 = 0.0;
931  this->a2 = 0.0;
932  enabled = false;
933  } else {
934  // Otherwise, increment the values
935  this->a1 += a1;
936  this->a2 += a2;
937  enabled = true;
938  }
939  this->reference = ref;
940  }
941 
942  bool is_enabled() const {
943  return enabled;
944  };
945 
946  void to_json(rapidjson::Value& el, rapidjson::Document& doc) {
947  el.AddMember("type", "IdealHelmholtzEnthalpyEntropyOffset", doc.GetAllocator());
948  el.AddMember("a1", static_cast<double>(a1), doc.GetAllocator());
949  el.AddMember("a2", static_cast<double>(a2), doc.GetAllocator());
950  };
951  void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
952 };
953 
959 class IdealHelmholtzLogTau : public BaseHelmholtzTerm
960 {
961  private:
962  CoolPropDbl a1;
963  bool enabled;
964 
965  public:
967  IdealHelmholtzLogTau() : a1(_HUGE), enabled(false) {}
968 
969  // Constructor
970  IdealHelmholtzLogTau(CoolPropDbl a1) : a1(a1), enabled(true) {}
971 
972  bool is_enabled() const {
973  return enabled;
974  };
975 
976  void to_json(rapidjson::Value& el, rapidjson::Document& doc) {
977  el.AddMember("type", "IdealHelmholtzLogTau", doc.GetAllocator());
978  el.AddMember("a1", static_cast<double>(a1), doc.GetAllocator());
979  };
980  void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
981 };
982 
988 class IdealHelmholtzPower : public BaseHelmholtzTerm
989 {
990 
991  private:
992  std::vector<CoolPropDbl> n, t; // Use these variables internally
993  std::size_t N;
994  bool enabled;
995 
996  public:
997  IdealHelmholtzPower() : N(0), enabled(false){};
998  // Constructor
999  IdealHelmholtzPower(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& t) : n(n), t(t), N(n.size()), enabled(true){};
1000 
1001  bool is_enabled() const {
1002  return enabled;
1003  };
1004 
1005  void to_json(rapidjson::Value& el, rapidjson::Document& doc) {
1006  el.AddMember("type", "IdealHelmholtzPower", doc.GetAllocator());
1007  cpjson::set_long_double_array("n", n, el, doc);
1008  cpjson::set_long_double_array("t", t, el, doc);
1009  };
1010  void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
1011 };
1012 
1059 class IdealHelmholtzPlanckEinsteinGeneralized : public BaseHelmholtzTerm
1060 {
1061 
1062  private:
1063  std::vector<CoolPropDbl> n, theta, c, d; // Use these variables internally
1064  std::size_t N;
1065  bool enabled;
1066 
1067  public:
1068  IdealHelmholtzPlanckEinsteinGeneralized() : N(0), enabled(false) {}
1069  // Constructor with std::vector instances
1070  IdealHelmholtzPlanckEinsteinGeneralized(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& theta,
1071  const std::vector<CoolPropDbl>& c, const std::vector<CoolPropDbl>& d)
1072  : n(n), theta(theta), c(c), d(d), N(n.size()), enabled(true) {}
1073 
1074  // Extend the vectors to allow for multiple instances feeding values to this function
1075  void extend(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& theta, const std::vector<CoolPropDbl>& c,
1076  const std::vector<CoolPropDbl>& d) {
1077  this->n.insert(this->n.end(), n.begin(), n.end());
1078  this->theta.insert(this->theta.end(), theta.begin(), theta.end());
1079  this->c.insert(this->c.end(), c.begin(), c.end());
1080  this->d.insert(this->d.end(), d.begin(), d.end());
1081  N += n.size();
1082  }
1083 
1084  bool is_enabled() const {
1085  return enabled;
1086  };
1087 
1088  void to_json(rapidjson::Value& el, rapidjson::Document& doc) {
1089  el.AddMember("type", "IdealHelmholtzPlanckEinsteinGeneralized", doc.GetAllocator());
1090  cpjson::set_long_double_array("n", n, el, doc);
1091  cpjson::set_long_double_array("theta", theta, el, doc);
1092  };
1093  void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
1094 };
1095 
1096 class IdealHelmholtzCP0Constant : public BaseHelmholtzTerm
1097 {
1098 
1099  private:
1100  double cp_over_R, Tc, T0, tau0; // Use these variables internally
1101  bool enabled;
1102 
1103  public:
1105  IdealHelmholtzCP0Constant() : cp_over_R(_HUGE), Tc(_HUGE), T0(_HUGE), tau0(_HUGE) {
1106  enabled = false;
1107  };
1108 
1110  IdealHelmholtzCP0Constant(CoolPropDbl cp_over_R, CoolPropDbl Tc, CoolPropDbl T0) : cp_over_R(cp_over_R), Tc(Tc), T0(T0) {
1111  enabled = true;
1112  tau0 = Tc / T0;
1113  };
1114 
1116  ~IdealHelmholtzCP0Constant(){};
1117 
1118  bool is_enabled() const {
1119  return enabled;
1120  };
1121 
1122  void to_json(rapidjson::Value& el, rapidjson::Document& doc) {
1123  el.AddMember("type", "IdealGasHelmholtzCP0Constant", doc.GetAllocator());
1124  el.AddMember("cp_over_R", cp_over_R, doc.GetAllocator());
1125  el.AddMember("Tc", Tc, doc.GetAllocator());
1126  el.AddMember("T0", T0, doc.GetAllocator());
1127  };
1128 
1129  void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
1130 };
1131 
1132 class IdealHelmholtzCP0PolyT : public BaseHelmholtzTerm
1133 {
1134  private:
1135  std::vector<CoolPropDbl> c, t;
1136  CoolPropDbl Tc, T0, tau0; // Use these variables internally
1137  std::size_t N;
1138  bool enabled;
1139 
1140  public:
1141  IdealHelmholtzCP0PolyT() : Tc(_HUGE), T0(_HUGE), tau0(_HUGE), N(0), enabled(false) {}
1142 
1144  IdealHelmholtzCP0PolyT(const std::vector<CoolPropDbl>& c, const std::vector<CoolPropDbl>& t, double Tc, double T0)
1145  : c(c), t(t), Tc(Tc), T0(T0), tau0(Tc / T0), N(c.size()), enabled(true) {
1146  assert(c.size() == t.size());
1147  }
1148 
1149  void extend(const std::vector<CoolPropDbl>& c, const std::vector<CoolPropDbl>& t) {
1150  this->c.insert(this->c.end(), c.begin(), c.end());
1151  this->t.insert(this->t.end(), t.begin(), t.end());
1152  N += c.size();
1153  }
1154 
1155  bool is_enabled() const {
1156  return enabled;
1157  };
1158 
1159  void to_json(rapidjson::Value& el, rapidjson::Document& doc);
1160  void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
1161 #if ENABLE_CATCH
1162  virtual mcx::MultiComplex<double> one_mcx(const mcx::MultiComplex<double>& tau, const mcx::MultiComplex<double>& delta) const override;
1163 #endif
1164 };
1168 class IdealHelmholtzGERG2004Sinh : public BaseHelmholtzTerm
1169 {
1170  private:
1171  std::vector<CoolPropDbl> n, theta;
1172  CoolPropDbl Tc, _Tr;
1173  std::size_t N;
1174  bool enabled;
1175 
1176  public:
1177  IdealHelmholtzGERG2004Sinh() : Tc(_HUGE), _Tr(_HUGE), N(0), enabled(false) {}
1178 
1180  IdealHelmholtzGERG2004Sinh(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& theta, double Tc)
1181  : n(n), theta(theta), Tc(Tc), _Tr(_HUGE), N(n.size()), enabled(true) {
1182  assert(n.size() == theta.size());
1183  }
1184 
1185  void extend(const std::vector<CoolPropDbl>& c, const std::vector<CoolPropDbl>& t) {
1186  this->n.insert(this->n.end(), n.begin(), n.end());
1187  this->theta.insert(this->theta.end(), theta.begin(), theta.end());
1188  N += c.size();
1189  }
1190  void set_Tred(CoolPropDbl Tr) {
1191  this->_Tr = Tr;
1192  }
1193 
1194  bool is_enabled() const {
1195  return enabled;
1196  };
1197  void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
1198 
1199  #if ENABLE_CATCH
1200  virtual mcx::MultiComplex<double> one_mcx(const mcx::MultiComplex<double>& tau, const mcx::MultiComplex<double>& delta) const override;
1201  #endif
1202 };
1203 
1204 class IdealHelmholtzGERG2004Cosh : public BaseHelmholtzTerm
1205 {
1206  private:
1207  std::vector<CoolPropDbl> n, theta;
1208  CoolPropDbl Tc, _Tr;
1209  std::size_t N;
1210  bool enabled;
1211 
1212  public:
1213  IdealHelmholtzGERG2004Cosh() : Tc(_HUGE), _Tr(_HUGE), N(0), enabled(false) {}
1214 
1216  IdealHelmholtzGERG2004Cosh(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& theta, double Tc)
1217  : n(n), theta(theta), Tc(Tc), _Tr(_HUGE), N(n.size()), enabled(true) {
1218  assert(n.size() == theta.size());
1219  }
1220 
1221  void extend(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& theta) {
1222  this->n.insert(this->n.end(), n.begin(), n.end());
1223  this->theta.insert(this->theta.end(), theta.begin(), theta.end());
1224  N += n.size();
1225  }
1226  void set_Tred(CoolPropDbl Tr) {
1227  this->_Tr = Tr;
1228  }
1229 
1230  bool is_enabled() const {
1231  return enabled;
1232  };
1233  void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
1234 
1235  #if ENABLE_CATCH
1236  virtual mcx::MultiComplex<double> one_mcx(const mcx::MultiComplex<double>& tau, const mcx::MultiComplex<double>& delta) const override;
1237  #endif
1238 
1239 };
1240 
1243 //\f[
1244 //\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
1245 //\f]
1246 //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
1247 //\f[
1248 //\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
1249 //\f]
1250 //or in terms of \f$ \tau \f$:
1251 //\f[
1252 //\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
1253 //\f]
1254 //Third partial:
1255 //\f[
1256 //\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)}
1257 //\f]
1258 //Now coming back to the ideal gas Helmholtz energy definition:
1259 //\f[
1260 //\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
1261 //\f]
1262 //Applying derivative
1263 //\f[
1264 //\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]
1265 //\f]
1266 //Fundamental theorem of calculus
1267 //\f[
1268 //\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}
1269 //\f]
1270 //Last two terms cancel, leaving
1271 //\f[
1272 //\alpha^0_{\tau} = -\int_{\tau_0}^{\tau} \frac{1}{\tau^2}\frac{c_p^0}{R_u}d\tau
1273 //\f]
1274 //Another derivative yields (from fundamental theorem of calculus)
1275 //\f[
1276 //\alpha^0_{\tau\tau} = - \frac{1}{\tau^2}\frac{c_p^0}{R_u}
1277 //\f]
1278 //
1279 //see also Jaeschke and Schley, 1995, (http://link.springer.com/article/10.1007%2FBF02083547#page-1)
1280 //*/
1282 //class IdealHelmholtzCP0AlyLee : public BaseHelmholtzTerm{
1283 //private:
1284 // std::vector<CoolPropDbl> c;
1285 // CoolPropDbl Tc, tau0, T0; // Use these variables internally
1286 // bool enabled;
1287 //public:
1288 // IdealHelmholtzCP0AlyLee(){enabled = false;};
1289 //
1290 // /// Constructor with std::vectors
1291 // IdealHelmholtzCP0AlyLee(const std::vector<CoolPropDbl> &c, double Tc, double T0)
1292 // :c(c), Tc(Tc), T0(T0)
1293 // {
1294 // tau0=Tc/T0;
1295 // enabled = true;
1296 // };
1297 //
1298 // /// Destructor
1299 // ~IdealHelmholtzCP0AlyLee(){};
1300 //
1301 // bool is_enabled() const {return enabled;};
1302 //
1303 // void to_json(rapidjson::Value &el, rapidjson::Document &doc);
1304 //
1305 //
1306 // /// The antiderivative given by \f$ \displaystyle\int \frac{1}{\tau^2}\frac{c_p^0}{R_u}d\tau \f$
1307 // /**
1308 // sympy code for this derivative:
1309 //
1310 // from sympy import *
1311 // a1,a2,a3,a4,a5,Tc,tau = symbols('a1,a2,a3,a4,a5,Tc,tau', real = True)
1312 // integrand = a1 + a2*(a3/Tc/sinh(a3*tau/Tc))**2 + a4*(a5/Tc/cosh(a5*tau/Tc))**2
1313 // integrand = integrand.rewrite(exp)
1314 // antideriv = trigsimp(integrate(integrand,tau))
1315 // display(antideriv)
1316 // print latex(antideriv)
1317 // print ccode(antideriv)
1318 //
1319 // \f[
1320 // \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]}
1321 // \f]
1322 // */
1323 // CoolPropDbl anti_deriv_cp0_tau2(const CoolPropDbl &tau);
1324 //
1325 // /// The antiderivative given by \f$ \displaystyle\int \frac{1}{\tau}\frac{c_p^0}{R_u}d\tau \f$
1326 // /**
1327 // sympy code for this derivative:
1328 //
1329 // 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)
1330 // 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
1331 //
1332 // term2 = a_1/tau*(a_2*tau/Tc/sinh(a_2*tau/Tc))**2
1333 // term2 = term2.rewrite(exp) # Unpack the sinh to exp functions
1334 // antideriv2 = trigsimp(integrate(term2,tau))
1335 // display(antideriv2)
1336 // print latex(antideriv2)
1337 // print ccode(antideriv2)
1338 //
1339 // term3 = a_3/tau*(a_4*tau/Tc/cosh(a_4*tau/Tc))**2
1340 // term3 = term3.rewrite(exp) # Unpack the cosh to exp functions
1341 // antideriv3 = factor(trigsimp(integrate(term3,tau).rewrite(exp)))
1342 // display(antideriv3)
1343 // print latex(antideriv3)
1344 // print ccode(antideriv3)
1345 //
1346 // Can be broken into three parts (trick is to express \f$sinh\f$ and \f$cosh\f$ in terms of \f$exp\f$ function)
1347 //
1348 // Term 2:
1349 // \f[
1350 // \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
1351 // \f]
1352 //
1353 // Term 3:
1354 // \f[
1355 // \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)
1356 // \f]
1357 // */
1358 // CoolPropDbl anti_deriv_cp0_tau(const CoolPropDbl &tau);
1359 //
1360 // CoolPropDbl base(const CoolPropDbl &tau, const CoolPropDbl &delta) throw();
1361 // CoolPropDbl dDelta(const CoolPropDbl &tau, const CoolPropDbl &delta) throw(){return 0.0;};
1362 // CoolPropDbl dTau(const CoolPropDbl &tau, const CoolPropDbl &delta) throw();
1363 // CoolPropDbl dDelta2(const CoolPropDbl &tau, const CoolPropDbl &delta) throw(){return 0.0;};
1364 // CoolPropDbl dDelta_dTau(const CoolPropDbl &tau, const CoolPropDbl &delta) throw(){return 0.0;};
1365 // CoolPropDbl dTau2(const CoolPropDbl &tau, const CoolPropDbl &delta) throw();
1366 // CoolPropDbl dDelta3(const CoolPropDbl &tau, const CoolPropDbl &delta) throw(){return 0.0;};
1367 // CoolPropDbl dDelta2_dTau(const CoolPropDbl &tau, const CoolPropDbl &delta) throw(){return 0.0;};
1368 // CoolPropDbl dDelta_dTau2(const CoolPropDbl &tau, const CoolPropDbl &delta) throw(){return 0.0;};
1369 // CoolPropDbl dTau3(const CoolPropDbl &tau, const CoolPropDbl &delta) throw();
1370 // CoolPropDbl dTau4(const CoolPropDbl &tau, const CoolPropDbl &delta) throw();
1371 //
1372 //};
1373 
1374 class IdealHelmholtzContainer : public BaseHelmholtzContainer
1375 {
1376  private:
1377  double _prefactor;
1378 
1379  public:
1380  IdealHelmholtzLead Lead;
1381  IdealHelmholtzEnthalpyEntropyOffset EnthalpyEntropyOffsetCore, EnthalpyEntropyOffset;
1382  IdealHelmholtzLogTau LogTau;
1383  IdealHelmholtzPower Power;
1384  IdealHelmholtzPlanckEinsteinGeneralized PlanckEinstein;
1385 
1386  IdealHelmholtzCP0Constant CP0Constant;
1387  IdealHelmholtzCP0PolyT CP0PolyT;
1388  IdealHelmholtzGERG2004Cosh GERG2004Cosh;
1389  IdealHelmholtzGERG2004Sinh GERG2004Sinh;
1390 
1391  IdealHelmholtzContainer() : _prefactor(1.0){};
1392 
1393  void set_prefactor(double prefactor) {
1394  _prefactor = prefactor;
1395  }
1396 
1397  void set_Tred(double T_red) {
1398  GERG2004Cosh.set_Tred(T_red);
1399  GERG2004Sinh.set_Tred(T_red);
1400  }
1401 
1402  void empty_the_EOS() {
1403  Lead = IdealHelmholtzLead();
1404  EnthalpyEntropyOffsetCore = IdealHelmholtzEnthalpyEntropyOffset();
1405  EnthalpyEntropyOffset = IdealHelmholtzEnthalpyEntropyOffset();
1406  LogTau = IdealHelmholtzLogTau();
1407  Power = IdealHelmholtzPower();
1408  PlanckEinstein = IdealHelmholtzPlanckEinsteinGeneralized();
1409  CP0Constant = IdealHelmholtzCP0Constant();
1410  CP0PolyT = IdealHelmholtzCP0PolyT();
1411  GERG2004Cosh = IdealHelmholtzGERG2004Cosh();
1412  GERG2004Sinh = IdealHelmholtzGERG2004Sinh();
1413  };
1414 
1415  HelmholtzDerivatives all(const CoolPropDbl tau, const CoolPropDbl delta, bool cache_values = false) {
1416  HelmholtzDerivatives derivs; // zeros out the elements
1417  Lead.all(tau, delta, derivs);
1418  EnthalpyEntropyOffsetCore.all(tau, delta, derivs);
1419  EnthalpyEntropyOffset.all(tau, delta, derivs);
1420  LogTau.all(tau, delta, derivs);
1421  Power.all(tau, delta, derivs);
1422  PlanckEinstein.all(tau, delta, derivs);
1423  CP0Constant.all(tau, delta, derivs);
1424  CP0PolyT.all(tau, delta, derivs);
1425  GERG2004Cosh.all(tau, delta, derivs);
1426  GERG2004Sinh.all(tau, delta, derivs);
1427 
1428  if (cache_values) {
1429  _base = derivs.alphar * _prefactor;
1430  _dDelta = derivs.dalphar_ddelta * _prefactor;
1431  _dTau = derivs.dalphar_dtau * _prefactor;
1432  _dDelta2 = derivs.d2alphar_ddelta2 * _prefactor;
1433  _dTau2 = derivs.d2alphar_dtau2 * _prefactor;
1434  _dDelta_dTau = derivs.d2alphar_ddelta_dtau * _prefactor;
1435  _dDelta3 = derivs.d3alphar_ddelta3 * _prefactor;
1436  _dTau3 = derivs.d3alphar_dtau3 * _prefactor;
1437  _dDelta2_dTau = derivs.d3alphar_ddelta2_dtau * _prefactor;
1438  _dDelta_dTau2 = derivs.d3alphar_ddelta_dtau2 * _prefactor;
1439  }
1440  return derivs * _prefactor;
1441  };
1442 };
1443 }; /* namespace CoolProp */
1444 
1445 #endif
Generated on Mon Mar 17 2025 04:25:23 for CoolProp by doxygen 1.9.1