1 #ifndef EXCESSHE_FUNCTIONS_H 2 #define EXCESSHE_FUNCTIONS_H 6 #include "CoolPropFluid.h" 7 #include "crossplatform_shared_ptr.h" 9 #include "Backends/Helmholtz/HelmholtzEOSMixtureBackend.h" 13 typedef std::vector<std::vector<CoolPropDbl>> STLMatrix;
33 virtual void update(
double tau,
double delta) {
35 phi.all(tau, delta, derivs);
37 double get(std::size_t itau, std::size_t idelta) {
38 return derivs.
get(itau, idelta);
43 phi.all(tau, delta, _derivs);
49 double dalphar_dDelta() {
50 return derivs.dalphar_ddelta;
52 double dalphar_dTau() {
53 return derivs.dalphar_dtau;
56 double d2alphar_dDelta2() {
57 return derivs.d2alphar_ddelta2;
59 double d2alphar_dDelta_dTau() {
60 return derivs.d2alphar_ddelta_dtau;
62 double d2alphar_dTau2() {
63 return derivs.d2alphar_dtau2;
66 double d3alphar_dTau3() {
67 return derivs.d3alphar_dtau3;
69 double d3alphar_dDelta_dTau2() {
70 return derivs.d3alphar_ddelta_dtau2;
72 double d3alphar_dDelta2_dTau() {
73 return derivs.d3alphar_ddelta2_dtau;
75 double d3alphar_dDelta3() {
76 return derivs.d3alphar_ddelta3;
79 double d4alphar_dTau4() {
80 return derivs.d4alphar_dtau4;
82 double d4alphar_dDelta_dTau3() {
83 return derivs.d4alphar_ddelta_dtau3;
85 double d4alphar_dDelta2_dTau2() {
86 return derivs.d4alphar_ddelta2_dtau2;
88 double d4alphar_dDelta3_dTau() {
89 return derivs.d4alphar_ddelta3_dtau;
91 double d4alphar_dDelta4() {
92 return derivs.d4alphar_ddelta4;
109 const std::vector<double>& eta,
const std::vector<double>& epsilon,
const std::vector<double>& beta,
110 const std::vector<double>& gamma, std::size_t Npower) {
113 std::vector<CoolPropDbl> _n(n.begin(), n.begin() + Npower);
114 std::vector<CoolPropDbl> _d(d.begin(), d.begin() + Npower);
115 std::vector<CoolPropDbl> _t(t.begin(), t.begin() + Npower);
116 std::vector<CoolPropDbl> _l(Npower, 0.0);
117 phi.add_Power(_n, _d, _t, _l);
119 if (n.size() == Npower) {
121 std::vector<CoolPropDbl> _n(n.begin() + Npower, n.end());
122 std::vector<CoolPropDbl> _d(d.begin() + Npower, d.end());
123 std::vector<CoolPropDbl> _t(t.begin() + Npower, t.end());
124 std::vector<CoolPropDbl> _eta(eta.begin() + Npower, eta.end());
125 std::vector<CoolPropDbl> _epsilon(epsilon.begin() + Npower, epsilon.end());
126 std::vector<CoolPropDbl> _beta(beta.begin() + Npower, beta.end());
127 std::vector<CoolPropDbl> _gamma(gamma.begin() + Npower, gamma.end());
128 phi.add_GERG2008Gaussian(_n, _d, _t, _eta, _epsilon, _beta, _gamma);
147 const std::vector<double>& l,
const std::vector<double>& eta,
const std::vector<double>& epsilon,
148 const std::vector<double>& beta,
const std::vector<double>& gamma, std::size_t Npower) {
151 std::vector<CoolPropDbl> _n(n.begin(), n.begin() + Npower);
152 std::vector<CoolPropDbl> _d(d.begin(), d.begin() + Npower);
153 std::vector<CoolPropDbl> _t(t.begin(), t.begin() + Npower);
154 std::vector<CoolPropDbl> _l(l.begin(), l.begin() + Npower);
155 phi.add_Power(_n, _d, _t, _l);
157 if (n.size() == Npower) {
159 std::vector<CoolPropDbl> _n(n.begin() + Npower, n.end());
160 std::vector<CoolPropDbl> _d(d.begin() + Npower, d.end());
161 std::vector<CoolPropDbl> _t(t.begin() + Npower, t.end());
162 std::vector<CoolPropDbl> _eta(eta.begin() + Npower, eta.end());
163 std::vector<CoolPropDbl> _epsilon(epsilon.begin() + Npower, epsilon.end());
164 std::vector<CoolPropDbl> _beta(beta.begin() + Npower, beta.end());
165 std::vector<CoolPropDbl> _gamma(gamma.begin() + Npower, gamma.end());
166 phi.add_Gaussian(_n, _d, _t, _eta, _epsilon, _beta, _gamma);
186 const std::vector<double>& l) {
187 std::vector<CoolPropDbl> _n(n.begin(), n.begin() + n.size());
188 std::vector<CoolPropDbl> _d(d.begin(), d.begin() + d.size());
189 std::vector<CoolPropDbl> _t(t.begin(), t.begin() + t.size());
190 std::vector<CoolPropDbl> _l(l.begin(), l.begin() + l.size());
191 phi.add_Power(_n, _d, _t, _l);
196 typedef shared_ptr<DepartureFunction> DepartureFunctionPointer;
202 std::vector<std::vector<DepartureFunctionPointer>> DepartureFunctionMatrix;
209 for (std::size_t i = 0; i < N; ++i) {
210 for (std::size_t j = 0; j < N; ++j) {
212 *(other.DepartureFunctionMatrix[i][j].get()) = *(other.DepartureFunctionMatrix[i][j].get());
222 for (std::size_t i = 0; i < N; ++i) {
223 for (std::size_t j = 0; j < N; ++j) {
225 _term.DepartureFunctionMatrix[i][j].reset(DepartureFunctionMatrix[i][j].
get()->copy_ptr());
236 F.resize(N, std::vector<CoolPropDbl>(N, 0));
237 DepartureFunctionMatrix.resize(N);
238 for (std::size_t i = 0; i < N; ++i) {
239 DepartureFunctionMatrix[i].resize(N);
244 for (std::size_t i = 0; i < N; i++) {
245 for (std::size_t j = i + 1; j < N; j++) {
246 DepartureFunctionMatrix[i][j]->update(tau, delta);
248 for (std::size_t j = 0; j < i; j++) {
249 DepartureFunctionMatrix[i][j]->update(tau, delta);
255 virtual HelmholtzDerivatives all(
const CoolPropDbl tau,
const CoolPropDbl delta,
const std::vector<CoolPropDbl>& mole_fractions,
256 bool cache_values =
false) {
264 if (cache_values ==
true) {
268 derivs.alphar = alphar(mole_fractions);
269 derivs.dalphar_ddelta = dalphar_dDelta(mole_fractions);
270 derivs.dalphar_dtau = dalphar_dTau(mole_fractions);
272 derivs.d2alphar_ddelta2 = d2alphar_dDelta2(mole_fractions);
273 derivs.d2alphar_ddelta_dtau = d2alphar_dDelta_dTau(mole_fractions);
274 derivs.d2alphar_dtau2 = d2alphar_dTau2(mole_fractions);
276 derivs.d3alphar_ddelta3 = d3alphar_dDelta3(mole_fractions);
277 derivs.d3alphar_ddelta2_dtau = d3alphar_dDelta2_dTau(mole_fractions);
278 derivs.d3alphar_ddelta_dtau2 = d3alphar_dDelta_dTau2(mole_fractions);
279 derivs.d3alphar_dtau3 = d3alphar_dTau3(mole_fractions);
281 derivs.d4alphar_ddelta4 = d4alphar_dDelta4(mole_fractions);
282 derivs.d4alphar_ddelta3_dtau = d4alphar_dDelta3_dTau(mole_fractions);
283 derivs.d4alphar_ddelta2_dtau2 = d4alphar_dDelta2_dTau2(mole_fractions);
284 derivs.d4alphar_ddelta_dtau3 = d4alphar_dDelta_dTau3(mole_fractions);
285 derivs.d4alphar_dtau4 = d4alphar_dTau4(mole_fractions);
288 return get_deriv_nocomp_notcached(mole_fractions, tau, delta);
291 HelmholtzDerivatives get_deriv_nocomp_notcached(
const std::vector<CoolPropDbl>& x,
double tau,
double delta)
const {
297 for (std::size_t i = 0; i < N - 1; i++) {
298 for (std::size_t j = i + 1; j < N; j++) {
300 DepartureFunctionMatrix[i][j]->calc_nocache(tau, delta, term);
301 summer = summer + term * x[i] * x[j] * F[i][j];
306 double get_deriv_nocomp_cached(
const std::vector<CoolPropDbl>& x, std::size_t itau, std::size_t idelta) {
312 for (std::size_t i = 0; i < N - 1; i++) {
313 for (std::size_t j = i + 1; j < N; j++) {
315 summer += x[i] * x[j] * F[i][j] * DepartureFunctionMatrix[i][j]->
get(itau, idelta);
320 double alphar(
const std::vector<CoolPropDbl>& x) {
321 return get_deriv_nocomp_cached(x, 0, 0);
323 double dalphar_dDelta(
const std::vector<CoolPropDbl>& x) {
324 return get_deriv_nocomp_cached(x, 0, 1);
326 double d2alphar_dDelta2(
const std::vector<CoolPropDbl>& x) {
327 return get_deriv_nocomp_cached(x, 0, 2);
329 double d2alphar_dDelta_dTau(
const std::vector<CoolPropDbl>& x) {
330 return get_deriv_nocomp_cached(x, 1, 1);
332 double dalphar_dTau(
const std::vector<CoolPropDbl>& x) {
333 return get_deriv_nocomp_cached(x, 1, 0);
335 double d2alphar_dTau2(
const std::vector<CoolPropDbl>& x) {
336 return get_deriv_nocomp_cached(x, 2, 0);
338 double d3alphar_dTau3(
const std::vector<CoolPropDbl>& x) {
339 return get_deriv_nocomp_cached(x, 3, 0);
341 double d3alphar_dDelta_dTau2(
const std::vector<CoolPropDbl>& x) {
342 return get_deriv_nocomp_cached(x, 2, 1);
344 double d3alphar_dDelta2_dTau(
const std::vector<CoolPropDbl>& x) {
345 return get_deriv_nocomp_cached(x, 1, 2);
347 double d3alphar_dDelta3(
const std::vector<CoolPropDbl>& x) {
348 return get_deriv_nocomp_cached(x, 0, 3);
350 double d4alphar_dTau4(
const std::vector<CoolPropDbl>& x) {
351 return get_deriv_nocomp_cached(x, 4, 0);
353 double d4alphar_dDelta_dTau3(
const std::vector<CoolPropDbl>& x) {
354 return get_deriv_nocomp_cached(x, 3, 1);
356 double d4alphar_dDelta2_dTau2(
const std::vector<CoolPropDbl>& x) {
357 return get_deriv_nocomp_cached(x, 2, 2);
359 double d4alphar_dDelta3_dTau(
const std::vector<CoolPropDbl>& x) {
360 return get_deriv_nocomp_cached(x, 1, 3);
362 double d4alphar_dDelta4(
const std::vector<CoolPropDbl>& x) {
363 return get_deriv_nocomp_cached(x, 0, 4);
366 double dalphar_dxi(
const std::vector<CoolPropDbl>& x, std::size_t i,
x_N_dependency_flag xN_flag) {
373 for (std::size_t k = 0; k < N; k++) {
375 summer += x[k] * F[i][k] * DepartureFunctionMatrix[i][k]->alphar();
383 CoolPropDbl dar_dxi = 0.0;
384 double FiNariN = F[i][N - 1] * DepartureFunctionMatrix[i][N - 1]->alphar();
385 dar_dxi += (1 - 2 * x[i]) * FiNariN;
386 for (std::size_t k = 0; k < N - 1; ++k) {
387 if (i == k)
continue;
388 double Fikarik = F[i][k] * DepartureFunctionMatrix[i][k]->alphar();
389 double FkNarkN = F[k][N - 1] * DepartureFunctionMatrix[k][N - 1]->alphar();
390 dar_dxi += x[k] * (Fikarik - FiNariN - FkNarkN);
394 throw ValueError(format(
"xN_flag is invalid"));
397 double d2alphardxidxj(
const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j,
x_N_dependency_flag xN_flag) {
404 return F[i][j] * DepartureFunctionMatrix[i][j]->alphar();
412 std::size_t N = x.size();
413 if (i == N - 1 || j == N - 1) {
416 double FiNariN = F[i][N - 1] * DepartureFunctionMatrix[i][N - 1]->alphar();
420 double Fijarij = F[i][j] * DepartureFunctionMatrix[i][j]->alphar();
421 double FjNarjN = F[j][N - 1] * DepartureFunctionMatrix[j][N - 1]->alphar();
422 return Fijarij - FiNariN - FjNarjN;
424 throw ValueError(format(
"xN_flag is invalid"));
427 double d3alphar_dxi_dxj_dDelta(
const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j,
x_N_dependency_flag xN_flag) {
434 return F[i][j] * DepartureFunctionMatrix[i][j]->dalphar_dDelta();
442 std::size_t N = x.size();
443 if (i == N - 1 || j == N - 1) {
446 double FiNariN = F[i][N - 1] * DepartureFunctionMatrix[i][N - 1]->dalphar_dDelta();
450 double Fijarij = F[i][j] * DepartureFunctionMatrix[i][j]->dalphar_dDelta();
451 double FjNarjN = F[j][N - 1] * DepartureFunctionMatrix[j][N - 1]->dalphar_dDelta();
452 return Fijarij - FiNariN - FjNarjN;
454 throw ValueError(format(
"xN_flag is invalid"));
457 double d3alphar_dxi_dxj_dTau(
const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j,
x_N_dependency_flag xN_flag) {
464 return F[i][j] * DepartureFunctionMatrix[i][j]->dalphar_dTau();
469 throw ValueError(format(
"xN_flag is invalid"));
472 double d4alphar_dxi_dxj_dDelta2(
const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j,
x_N_dependency_flag xN_flag) {
479 return F[i][j] * DepartureFunctionMatrix[i][j]->d2alphar_dDelta2();
484 throw ValueError(format(
"xN_flag is invalid"));
487 double d4alphar_dxi_dxj_dDelta_dTau(
const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j,
x_N_dependency_flag xN_flag) {
494 return F[i][j] * DepartureFunctionMatrix[i][j]->d2alphar_dDelta_dTau();
502 double FiNariN = F[i][N - 1] * DepartureFunctionMatrix[i][N - 1]->d2alphar_dDelta_dTau();
503 CoolPropDbl d3ar_dxi_dDelta_dTau = (1 - 2 * x[i]) * FiNariN;
504 for (std::size_t k = 0; k < N - 1; ++k) {
505 if (i == k)
continue;
506 double Fikarik = F[i][k] * DepartureFunctionMatrix[i][k]->d2alphar_dDelta_dTau();
507 double FkNarkN = F[k][N - 1] * DepartureFunctionMatrix[k][N - 1]->d2alphar_dDelta_dTau();
508 d3ar_dxi_dDelta_dTau += x[k] * (Fikarik - FiNariN - FkNarkN);
510 return d3ar_dxi_dDelta_dTau;
512 throw ValueError(format(
"xN_flag is invalid"));
515 double d4alphar_dxi_dxj_dTau2(
const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j,
x_N_dependency_flag xN_flag) {
522 return F[i][j] * DepartureFunctionMatrix[i][j]->d2alphar_dTau2();
527 throw ValueError(format(
"xN_flag is invalid"));
531 double d3alphardxidxjdxk(
const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j, std::size_t k,
x_N_dependency_flag xN_flag) {
534 double d2alphar_dxi_dTau(
const std::vector<CoolPropDbl>& x, std::size_t i,
x_N_dependency_flag xN_flag) {
541 for (std::size_t k = 0; k < N; k++) {
543 summer += x[k] * F[i][k] * DepartureFunctionMatrix[i][k]->dalphar_dTau();
551 double FiNariN = F[i][N - 1] * DepartureFunctionMatrix[i][N - 1]->dalphar_dTau();
552 CoolPropDbl d2ar_dxi_dTau = (1 - 2 * x[i]) * FiNariN;
553 for (std::size_t k = 0; k < N - 1; ++k) {
554 if (i == k)
continue;
555 double Fikarik = F[i][k] * DepartureFunctionMatrix[i][k]->dalphar_dTau();
556 double FkNarkN = F[k][N - 1] * DepartureFunctionMatrix[k][N - 1]->dalphar_dTau();
557 d2ar_dxi_dTau += x[k] * (Fikarik - FiNariN - FkNarkN);
559 return d2ar_dxi_dTau;
561 throw ValueError(format(
"xN_flag is invalid"));
564 double d2alphar_dxi_dDelta(
const std::vector<CoolPropDbl>& x, std::size_t i,
x_N_dependency_flag xN_flag) {
571 for (std::size_t k = 0; k < N; k++) {
573 summer += x[k] * F[i][k] * DepartureFunctionMatrix[i][k]->dalphar_dDelta();
581 CoolPropDbl d2ar_dxi_dDelta = 0;
582 double FiNariN = F[i][N - 1] * DepartureFunctionMatrix[i][N - 1]->dalphar_dDelta();
583 d2ar_dxi_dDelta += (1 - 2 * x[i]) * FiNariN;
584 for (std::size_t k = 0; k < N - 1; ++k) {
585 if (i == k)
continue;
586 double Fikarik = F[i][k] * DepartureFunctionMatrix[i][k]->dalphar_dDelta();
587 double FkNarkN = F[k][N - 1] * DepartureFunctionMatrix[k][N - 1]->dalphar_dDelta();
588 d2ar_dxi_dDelta += x[k] * (Fikarik - FiNariN - FkNarkN);
590 return d2ar_dxi_dDelta;
592 throw ValueError(format(
"xN_flag is invalid"));
595 double d3alphar_dxi_dDelta2(
const std::vector<CoolPropDbl>& x, std::size_t i,
x_N_dependency_flag xN_flag) {
602 for (std::size_t k = 0; k < N; k++) {
604 summer += x[k] * F[i][k] * DepartureFunctionMatrix[i][k]->d2alphar_dDelta2();
612 double FiNariN = F[i][N - 1] * DepartureFunctionMatrix[i][N - 1]->d2alphar_dDelta2();
613 CoolPropDbl d3ar_dxi_dDelta2 = (1 - 2 * x[i]) * FiNariN;
614 for (std::size_t k = 0; k < N - 1; ++k) {
615 if (i == k)
continue;
616 double Fikarik = F[i][k] * DepartureFunctionMatrix[i][k]->d2alphar_dDelta2();
617 double FkNarkN = F[k][N - 1] * DepartureFunctionMatrix[k][N - 1]->d2alphar_dDelta2();
618 d3ar_dxi_dDelta2 += x[k] * (Fikarik - FiNariN - FkNarkN);
620 return d3ar_dxi_dDelta2;
622 throw ValueError(format(
"xN_flag is invalid"));
625 double d4alphar_dxi_dDelta3(
const std::vector<CoolPropDbl>& x, std::size_t i,
x_N_dependency_flag xN_flag) {
632 for (std::size_t k = 0; k < N; k++) {
634 summer += x[k] * F[i][k] * DepartureFunctionMatrix[i][k]->d3alphar_dDelta3();
639 throw ValueError(format(
"xN_flag is invalid"));
642 double d3alphar_dxi_dTau2(
const std::vector<CoolPropDbl>& x, std::size_t i,
x_N_dependency_flag xN_flag) {
649 for (std::size_t k = 0; k < N; k++) {
651 summer += x[k] * F[i][k] * DepartureFunctionMatrix[i][k]->d2alphar_dTau2();
659 double FiNariN = F[i][N - 1] * DepartureFunctionMatrix[i][N - 1]->d2alphar_dTau2();
660 CoolPropDbl d3ar_dxi_dTau2 = (1 - 2 * x[i]) * FiNariN;
661 for (std::size_t k = 0; k < N - 1; ++k) {
662 if (i == k)
continue;
663 double Fikarik = F[i][k] * DepartureFunctionMatrix[i][k]->d2alphar_dTau2();
664 double FkNarkN = F[k][N - 1] * DepartureFunctionMatrix[k][N - 1]->d2alphar_dTau2();
665 d3ar_dxi_dTau2 += x[k] * (Fikarik - FiNariN - FkNarkN);
667 return d3ar_dxi_dTau2;
669 throw ValueError(format(
"xN_flag is invalid"));
672 double d4alphar_dxi_dTau3(
const std::vector<CoolPropDbl>& x, std::size_t i,
x_N_dependency_flag xN_flag) {
679 for (std::size_t k = 0; k < N; k++) {
681 summer += x[k] * F[i][k] * DepartureFunctionMatrix[i][k]->d3alphar_dTau3();
686 throw ValueError(format(
"xN_flag is invalid"));
689 double d3alphar_dxi_dDelta_dTau(
const std::vector<CoolPropDbl>& x, std::size_t i,
x_N_dependency_flag xN_flag) {
696 for (std::size_t k = 0; k < N; k++) {
698 summer += x[k] * F[i][k] * DepartureFunctionMatrix[i][k]->d2alphar_dDelta_dTau();
703 throw ValueError(format(
"xN_flag is invalid"));
706 double d4alphar_dxi_dDelta2_dTau(
const std::vector<CoolPropDbl>& x, std::size_t i,
x_N_dependency_flag xN_flag) {
713 for (std::size_t k = 0; k < N; k++) {
715 summer += x[k] * F[i][k] * DepartureFunctionMatrix[i][k]->d3alphar_dDelta2_dTau();
720 throw ValueError(format(
"xN_flag is invalid"));
723 double d4alphar_dxi_dDelta_dTau2(
const std::vector<CoolPropDbl>& x, std::size_t i,
x_N_dependency_flag xN_flag) {
730 for (std::size_t k = 0; k < N; k++) {
732 summer += x[k] * F[i][k] * DepartureFunctionMatrix[i][k]->d3alphar_dDelta_dTau2();
737 throw ValueError(format(
"xN_flag is invalid"));
A generalized residual helmholtz energy container that can deal with a wide range of terms which can ...
Definition: Helmholtz.h:332
double get(std::size_t itau, std::size_t idelta)
Retrieve a single value based on the number of derivatives with respect to tau and delta...
Definition: Helmholtz.h:77
A polynomial/exponential departure function.
Definition: ExcessHEFunction.h:181
GaussianExponentialDepartureFunction(const std::vector< double > &n, const std::vector< double > &d, const std::vector< double > &t, const std::vector< double > &l, const std::vector< double > &eta, const std::vector< double > &epsilon, const std::vector< double > &beta, const std::vector< double > &gamma, std::size_t Npower)
Definition: ExcessHEFunction.h:146
x_N is an dependent variable, calculated by
Definition: ReducingFunctions.h:21
void update(double tau, double delta)
Update the internal cached derivatives in each departure function.
Definition: ExcessHEFunction.h:243
Definition: ExcessHEFunction.h:198
Definition: Exceptions.h:45
void resize(std::size_t N)
Resize the parts of this term.
Definition: ExcessHEFunction.h:234
virtual HelmholtzDerivatives all(const CoolPropDbl tau, const CoolPropDbl delta, const std::vector< CoolPropDbl > &mole_fractions, bool cache_values=false)
Calculate all the derivatives that do not involve any composition derivatives.
Definition: ExcessHEFunction.h:255
GERG2008DepartureFunction(const std::vector< double > &n, const std::vector< double > &d, const std::vector< double > &t, const std::vector< double > &eta, const std::vector< double > &epsilon, const std::vector< double > &beta, const std::vector< double > &gamma, std::size_t Npower)
Definition: ExcessHEFunction.h:108
Definition: Helmholtz.h:45
x_N is an independent variable, and not calculated by
Definition: ReducingFunctions.h:20
A hybrid gaussian with temperature and density dependence along with.
Definition: ExcessHEFunction.h:142
x_N_dependency_flag
Definition: ReducingFunctions.h:18
This file contains flash routines in which the state is unknown, and a solver of some kind must be us...
Definition: AbstractState.h:19
The abstract base class for departure functions used in the excess part of the Helmholtz energy...
Definition: ExcessHEFunction.h:20
The departure function used by the GERG-2008 formulation.
Definition: ExcessHEFunction.h:104