From e5f2b31f3f80407bf4b272b2d3769a75f077103a Mon Sep 17 00:00:00 2001 From: David Andrs Date: Tue, 23 May 2023 10:13:44 -0600 Subject: [PATCH 1/5] Constant in the collision integral is a parameter --- include/TransportModels.h | 13 ++++++++++--- src/Air.cpp | 2 +- src/Nitrogen.cpp | 2 +- test/src/TransportModels_test.cpp | 2 +- 4 files changed, 13 insertions(+), 6 deletions(-) diff --git a/include/TransportModels.h b/include/TransportModels.h index d0ea913..e2b49d9 100644 --- a/include/TransportModels.h +++ b/include/TransportModels.h @@ -49,7 +49,7 @@ class Eta0AndPoly { /// /// @tparam T The basic data type /// -/// \f$ \eta_0(T) = \frac{0.0266958 \sqrt{M T}}{\sigma^2 \Omega(T^*)} \f$ +/// \f$ \eta_0(T) = \frac{C \sqrt{M T}}{\sigma^2 \Omega(T^*)} \f$ /// /// where \f$\sigma\f$ is the Lennard-Jones size parameter and /// \f$\Omega\f$ is the collision integral, given by @@ -63,11 +63,17 @@ class LennardJones { public: /// Lennard-Jones model /// + /// @param C Constant in front of term /// @param M Molar mass \f$[{kg\over mol}]\f$ /// @param epsilon_over_k \f${\epsilon\over k} [K]\f$ /// @param sigma \f$\sigma\f$ /// @param b \f$b_i\f$ - LennardJones(double M, double epsilon_over_k, double sigma, const std::vector & b) : + LennardJones(double C, + double M, + double epsilon_over_k, + double sigma, + const std::vector & b) : + C(C), M(M), epsilon_over_k(epsilon_over_k), sigma(sigma), @@ -87,11 +93,12 @@ class LennardJones { for (unsigned int i = 0; i < this->b.size(); i++) Omega_T_star += this->b[i] * std::pow(log_T_star, i); Omega_T_star = std::exp(Omega_T_star); - return 0.0266958 * std::sqrt(1000.0 * this->M * temperature) / + return this->C * std::sqrt(1000.0 * this->M * temperature) / (this->sigma * this->sigma * Omega_T_star); } private: + double C; double M; double epsilon_over_k; const double sigma; diff --git a/src/Air.cpp b/src/Air.cpp index 119c5eb..c1d8872 100644 --- a/src/Air.cpp +++ b/src/Air.cpp @@ -46,7 +46,7 @@ Air::Air() : { 1, 3, 5, 6, 1, 3, 11, 1, 3 }, { 1.6, 0.8, 0.95, 1.25, 3.6, 6, 3.25, 3.5, 15 }, { 1, 1, 1, 1, 2, 2, 2, 3, 3 }), - eta_0(MOLAR_MASS, 103.3, 0.360, { 0.431, -0.4623, 0.08406, 0.005341, -0.00331 }), + eta_0(0.0266958, MOLAR_MASS, 103.3, 0.360, { 0.431, -0.4623, 0.08406, 0.005341, -0.00331 }), eta_r({ 10.72, 1.122, 0.002019, -8.876, -0.02916 }, { 0.2, 0.05, 2.4, 0.6, 3.6 }, { 1, 4, 9, 1, 8 }, diff --git a/src/Nitrogen.cpp b/src/Nitrogen.cpp index 2d535e2..e72d4a9 100644 --- a/src/Nitrogen.cpp +++ b/src/Nitrogen.cpp @@ -38,7 +38,7 @@ Nitrogen::Nitrogen() : { 325, 325, 300, 275 }, { 1.16, 1.16, 1.13, 1.25 }), - eta_0(MOLAR_MASS, 98.94, 0.3656, { 0.431, -0.4623, 0.08406, 0.005341, -0.00331 }), + eta_0(0.0266958, MOLAR_MASS, 98.94, 0.3656, { 0.431, -0.4623, 0.08406, 0.005341, -0.00331 }), eta_r({ 10.72, 0.03989, 0.001208, -7.402, 4.62 }, { 0.1, 0.25, 3.2, 0.9, 0.3 }, { 2, 10, 12, 2, 1 }, diff --git a/test/src/TransportModels_test.cpp b/test/src/TransportModels_test.cpp index 222f6cd..196f5d2 100644 --- a/test/src/TransportModels_test.cpp +++ b/test/src/TransportModels_test.cpp @@ -11,7 +11,7 @@ TEST(TransportTest, eta0_and_poly) TEST(TransportTest, lennard_jones) { - LennardJones lj(2.e-3, 3., 2., { 2, 3 }); + LennardJones lj(0.0266958, 2.e-3, 3., 2., { 2, 3 }); EXPECT_DOUBLE_EQ(lj.value(2.), 0.0060967411665096916); } From 828202b1369d1ee90b10324c984666b7c56992d9 Mon Sep 17 00:00:00 2001 From: David Andrs Date: Mon, 22 May 2023 15:04:05 -0600 Subject: [PATCH 2/5] Adding non-analytic model to Helmholtz --- include/Helmholtz.h | 288 ++++++++++++++++++++++++++++++++++++ test/src/Helmholtz_test.cpp | 22 +++ 2 files changed, 310 insertions(+) diff --git a/include/Helmholtz.h b/include/Helmholtz.h index fdb5970..a204323 100644 --- a/include/Helmholtz.h +++ b/include/Helmholtz.h @@ -691,6 +691,294 @@ class Helmholtz : public SinglePhaseFluidProperties { std::vector beta; std::vector gamma; }; + + template + class ResidualNonAnalytic { + public: + ResidualNonAnalytic(const std::vector & n, + const std::vector & a, + const std::vector & b, + const std::vector & beta, + const std::vector & A, + const std::vector & B, + const std::vector & C, + const std::vector & D) : + n(n), + a(a), + b(b), + beta(beta), + A(A), + B(B), + C(C), + D(D) + { + assert(n.size() == a.size()); + assert(n.size() == b.size()); + assert(n.size() == beta.size()); + assert(n.size() == A.size()); + assert(n.size() == B.size()); + assert(n.size() == C.size()); + assert(n.size() == D.size()); + } + + T + alpha(T delta, T tau) const + { + fiddle(delta, tau); + T sum = 0; + for (std::size_t i = 0; i < this->n.size(); i++) + sum += delta * this->n[i] * DELTA_bi(i, delta, tau) * PSI(i, delta, tau); + return sum; + } + + T + ddelta(T delta, T tau) const + { + fiddle(delta, tau); + T sum = 0; + for (std::size_t i = 0; i < this->n.size(); i++) { + sum += this->n[i] * (DELTA_bi(i, delta, tau) * PSI(i, delta, tau) + + delta * dDELTA_bi_ddelta(i, delta, tau) * PSI(i, delta, tau) + + delta * DELTA_bi(i, delta, tau) * dPSI_ddelta(i, delta, tau)); + } + return sum; + } + + T + dtau(T delta, T tau) const + { + fiddle(delta, tau); + T sum = 0; + for (std::size_t i = 0; i < this->n.size(); i++) { + sum += this->n[i] * delta * + (dDELTA_bi_dtau(i, delta, tau) * PSI(i, delta, tau) + + DELTA_bi(i, delta, tau) * dPSI_dtau(i, delta, tau)); + } + return sum; + } + + T + d2delta(T delta, T tau) const + { + fiddle(delta, tau); + T sum = 0; + for (std::size_t i = 0; i < this->n.size(); i++) { + sum += this->n[i] * + (dDELTA_bi_ddelta(i, delta, tau) * PSI(i, delta, tau) + + DELTA_bi(i, delta, tau) * dPSI_ddelta(i, delta, tau) + + dDELTA_bi_ddelta(i, delta, tau) * PSI(i, delta, tau) + + delta * d2DELTA_bi_ddelta2(i, delta, tau) * PSI(i, delta, tau) + + delta * dDELTA_bi_ddelta(i, delta, tau) * dPSI_ddelta(i, delta, tau) + + DELTA_bi(i, delta, tau) * dPSI_ddelta(i, delta, tau) + + delta * dDELTA_bi_ddelta(i, delta, tau) * dPSI_ddelta(i, delta, tau) + + delta * DELTA_bi(i, delta, tau) * d2PSI_ddelta2(i, delta, tau)); + } + return sum; + } + + T + d2tau(T delta, T tau) const + { + fiddle(delta, tau); + T sum = 0; + for (std::size_t i = 0; i < this->n.size(); i++) { + sum += this->n[i] * delta * + (d2DELTA_bi_dtau2(i, delta, tau) * PSI(i, delta, tau) + + 2 * dDELTA_bi_dtau(i, delta, tau) * dPSI_dtau(i, delta, tau) + + DELTA_bi(i, delta, tau) * d2PSI_dtau2(i, delta, tau)); + } + return sum; + } + + T + d2deltatau(T delta, T tau) const + { + fiddle(delta, tau); + T sum = 0; + for (std::size_t i = 0; i < this->n.size(); i++) { + sum += this->n[i] * + (dDELTA_bi_dtau(i, delta, tau) * PSI(i, delta, tau) + + delta * (d2DELTA_bi_ddeltatau(i, delta, tau) * PSI(i, delta, tau) + + dDELTA_bi_dtau(i, delta, tau) * dPSI_ddelta(i, delta, tau) + + dDELTA_bi_ddelta(i, delta, tau) * dPSI_dtau(i, delta, tau) + + DELTA_bi(i, delta, tau) * d2PSI_ddeltatau(i, delta, tau))); + } + return sum; + } + + private: + // Manipulate `delta` and `tau` such that we are not evaluating at the critical point but + // very close to it + void + fiddle(T & delta, T & tau) const + { + T EPSILON = std::numeric_limits::epsilon(); + if (std::abs(tau - 1) < 10 * EPSILON) + tau = 1.0 + 10 * EPSILON; + if (std::abs(delta - 1) < 10 * EPSILON) + delta = 1.0 + 10 * EPSILON; + } + + T + theta(std::size_t i, T delta, T tau) const + { + return (1.0 - tau) + + this->A[i] * std::pow(sqr(delta - 1.0), 1.0 / (2.0 * this->beta[i])); + } + + T + dtheta_ddelta(std::size_t i, T delta, T tau) const + { + return this->A[i] * std::pow(sqr(delta - 1), 1. / (2. * this->beta[i])) / + (this->beta[i] * (delta - 1)); + } + + T + dtheta_dtau(std::size_t i, T delta, T tau) const + { + return -1; + } + + T + d2theta_ddelta2(std::size_t i, T delta, T tau) const + { + return -1 * + (this->A[i] * (this->beta[i] - 1) * + std::pow(sqr(delta - 1), 1. / (2 * this->beta[i]) - 1.)) / + (sqr(this->beta[i])); + } + + T + DELTA(std::size_t i, T delta, T tau) const + { + return sqr(theta(i, delta, tau)) + this->B[i] * std::pow(sqr(delta - 1.0), this->a[i]); + } + + T + dDELTA_ddelta(std::size_t i, T delta, T tau) const + { + return 2. * theta(i, delta, tau) * dtheta_ddelta(i, delta, tau) + + (2 * this->a[i] * this->B[i] * std::pow(delta - 1, 2 * this->a[i] - 1)); + } + + T + dDELTA_dtau(std::size_t i, T delta, T tau) const + { + return 2. * theta(i, delta, tau) * dtheta_dtau(i, delta, tau); + } + + T + d2DELTA_ddelta2(std::size_t i, T delta, T tau) const + { + return 2. * (sqr(dtheta_ddelta(i, delta, tau) + + theta(i, delta, tau) * d2theta_ddelta2(i, delta, tau))) + + 2 * this->a[i] * (2 * this->a[i] - 1) * this->B[i] * + std::pow(delta - 1, 2 * this->a[i] - 2); + } + + T + d2DELTA_dtau2(std::size_t i, T delta, T tau) const + { + return 2. * sqr(dtheta_dtau(i, delta, tau)); + } + + T + d2DELTA_ddeltatau(std::size_t i, T delta, T tau) const + { + return 2 * dtheta_ddelta(i, delta, tau) * dtheta_dtau(i, delta, tau); + } + + T + DELTA_bi(std::size_t i, T delta, T tau) const + { + return std::pow(DELTA(i, delta, tau), this->b[i]); + } + + T + dDELTA_bi_ddelta(std::size_t i, T delta, T tau) const + { + return this->b[i] * std::pow(DELTA(i, delta, tau), this->b[i] - 1) * + dDELTA_ddelta(i, delta, tau); + } + + T + dDELTA_bi_dtau(std::size_t i, T delta, T tau) const + { + return this->b[i] * std::pow(DELTA(i, delta, tau), this->b[i] - 1) * + dDELTA_dtau(i, delta, tau); + } + + T + d2DELTA_bi_ddelta2(std::size_t i, T delta, T tau) const + { + return this->b[i] * std::pow(DELTA(i, delta, tau), this->b[i] - 1) * + d2DELTA_ddelta2(i, delta, tau); + } + + T + d2DELTA_bi_dtau2(std::size_t i, T delta, T tau) const + { + return this->b[i] * std::pow(DELTA(i, delta, tau), this->b[i] - 1) * + d2DELTA_dtau2(i, delta, tau); + } + + T + d2DELTA_bi_ddeltatau(std::size_t i, T delta, T tau) const + { + return this->b[i] * std::pow(DELTA(i, delta, tau), this->b[i] - 1) * + d2DELTA_ddeltatau(i, delta, tau); + } + + T + PSI(std::size_t i, T delta, T tau) const + { + return std::exp(-this->C[i] * sqr(delta - 1.0) - this->D[i] * sqr(tau - 1.0)); + } + + T + dPSI_ddelta(std::size_t i, T delta, T tau) const + { + return -2 * this->C[i] * (delta - 1) * + std::exp(-this->C[i] * sqr(delta - 1) - this->D[i] * sqr(tau - 1)); + } + + T + dPSI_dtau(std::size_t i, T delta, T tau) const + { + return -2 * this->D[i] * (tau - 1) * + std::exp(-this->C[i] * sqr(delta - 1) - this->D[i] * sqr(tau - 1)); + } + + T + d2PSI_ddelta2(std::size_t i, T delta, T tau) const + { + return 2 * this->C[i] * (2 * this->C[i] * sqr(delta - 1) - 1) * + std::exp(-this->C[i] * sqr(delta - 1) - this->D[i] * sqr(tau - 1)); + } + + T + d2PSI_dtau2(std::size_t i, T delta, T tau) const + { + return 2 * this->D[i] * (2 * this->D[i] * sqr(tau - 1) - 1) * + std::exp(-this->C[i] * sqr(delta - 1) - this->D[i] * sqr(tau - 1)); + } + + T + d2PSI_ddeltatau(std::size_t i, T delta, T tau) const + { + return 4. * this->C[i] * (delta - 1) * this->D[i] * (tau - 1) * + std::exp(-this->C[i] * sqr(delta - 1) - this->D[i] * sqr(tau - 1)); + } + + std::vector n; + std::vector a; + std::vector b; + std::vector beta; + std::vector A; + std::vector B; + std::vector C; + std::vector D; + }; }; } // namespace fprops diff --git a/test/src/Helmholtz_test.cpp b/test/src/Helmholtz_test.cpp index f7601c1..299bb84 100644 --- a/test/src/Helmholtz_test.cpp +++ b/test/src/Helmholtz_test.cpp @@ -219,3 +219,25 @@ TEST(HelmholtzTest, residual_gaussian) Test t; } + +TEST(HelmholtzTest, non_analytic) +{ + class Test : public MockHelmholtz { + public: + Test() : + MockHelmholtz(), + a({ 2, 3 }, { 3, 4 }, { 4, 5 }, { 5, 6 }, { 3, 2 }, { 4, 3 }, { 5, 4 }, { 6, 5 }) + { + EXPECT_DOUBLE_EQ(a.alpha(2., 3.), 5.5677378067433475e-08); + EXPECT_DOUBLE_EQ(a.ddelta(2., 3.), 1.7956263727569099e-06); + EXPECT_DOUBLE_EQ(a.dtau(2., 3.), -1.1171086932549971e-06); + EXPECT_DOUBLE_EQ(a.d2delta(2., 3.), -1.5082537326175269e-05); + EXPECT_DOUBLE_EQ(a.d2tau(2., 3.), 2.2058154287500501e-05); + EXPECT_DOUBLE_EQ(a.d2deltatau(2., 3.), -3.5433699047353868e-05); + } + + Helmholtz::ResidualNonAnalytic a; + }; + + Test t; +} From f1b40b4d3842fe428e4a174ac2eb9221f427a849 Mon Sep 17 00:00:00 2001 From: David Andrs Date: Mon, 22 May 2023 15:05:48 -0600 Subject: [PATCH 3/5] Adding carbon dioxide fluid properties (c++) --- include/CarbonDioxide.h | 50 +++++++++ src/CMakeLists.txt | 1 + src/CarbonDioxide.cpp | 192 ++++++++++++++++++++++++++++++++ test/src/CMakeLists.txt | 1 + test/src/CarbonDioxide_test.cpp | 48 ++++++++ 5 files changed, 292 insertions(+) create mode 100644 include/CarbonDioxide.h create mode 100644 src/CarbonDioxide.cpp create mode 100644 test/src/CarbonDioxide_test.cpp diff --git a/include/CarbonDioxide.h b/include/CarbonDioxide.h new file mode 100644 index 0000000..f39c0c6 --- /dev/null +++ b/include/CarbonDioxide.h @@ -0,0 +1,50 @@ +#pragma once + +#include "Helmholtz.h" +#include "TransportModels.h" + +namespace fprops { + +/// Carbon dioxide fluid properties +/// +/// References: +/// 1. R. Span and W. Wagner. A New Equation of State for Carbon Dioxide Covering the Fluid Region +/// from the Triple Point Temperature to 1100 K at Pressures up to 800 MPa. J. Phys. Chem. Ref. +/// Data, 25:1509–1596, 1996. doi:10.1063/1.555991. +/// 2. G. Scalabrin, P. Marchi, F. Finezzo, and R. Span. A Reference Multiparameter Thermal +/// Conductivity Equation for Carbon Dioxide with an Optimized Functional Form. J. Phys. Chem. +/// Ref. Data, 35(4):1549–1575, 2006. doi:10.1063/1.2213631. +/// 3. A. Fenghour, W.A. Wakeham, and V. Vesovic. The viscosity of carbon dioxide. J. Phys. Chem. +/// Ref. Data, 27(1):31–44, 1998. 5. doi:10.1063/1.556013. +class CarbonDioxide : public Helmholtz { +public: + CarbonDioxide(); + +protected: + [[nodiscard]] double alpha(double delta, double tau) const override; + [[nodiscard]] double dalpha_ddelta(double delta, double tau) const override; + [[nodiscard]] double dalpha_dtau(double delta, double tau) const override; + [[nodiscard]] double d2alpha_ddelta2(double delta, double tau) const override; + [[nodiscard]] double d2alpha_dtau2(double delta, double tau) const override; + [[nodiscard]] double d2alpha_ddeltatau(double delta, double tau) const override; + + [[nodiscard]] double mu_from_rho_T(double rho, double T) const override; + [[nodiscard]] double k_from_rho_T(double rho, double T) const override; + +private: + IdealGasLead lead; + IdealGasLogTau log_tau; + IdealGasPlanckEinstein pe; + IdealEnthalpyEntropyOffset offset; + + ResidualPower power_r; + ResidualPowerExp power_exp_r; + ResidualGaussian gauss; + ResidualNonAnalytic noan; + + LennardJones eta_0; + ModifiedBatshinskiHildebrand eta_r; + ModifiedBatshinskiHildebrand lambda_r; +}; + +} // namespace fprops diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index d86dfe7..00b0d94 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -2,6 +2,7 @@ add_library( ${PROJECT_NAME} SHARED Air.cpp + CarbonDioxide.cpp FluidProperties.cpp Helium.cpp Helmholtz.cpp diff --git a/src/CarbonDioxide.cpp b/src/CarbonDioxide.cpp new file mode 100644 index 0000000..ed5ed7a --- /dev/null +++ b/src/CarbonDioxide.cpp @@ -0,0 +1,192 @@ +#include "CarbonDioxide.h" + +namespace fprops { + +static const double GAS_CONSTANT = 8.31451; +static const double MOLAR_MASS = 44.0098e-3; +static const double T_CRIT = 304.1282; +static const double RHO_MOLAR_CRIT = 10624.9063; +static const double RHO_CRIT = RHO_MOLAR_CRIT * MOLAR_MASS; + +CarbonDioxide::CarbonDioxide() : + Helmholtz(GAS_CONSTANT, MOLAR_MASS, RHO_CRIT, T_CRIT), + lead(8.37304456, -3.70454304), + log_tau(2.5), + pe({ 1.99427042, 0.62105248, 0.41195293, 1.04028922, 0.08327678 }, + { 3.15163, 6.1119, 6.77708, 11.32384, 27.08792 }), + offset(-14.4979156224319, 8.82013935801453), + power_r({ 0.388568232032, + 2.93854759427, + -5.5867188535, + -0.767531995925, + 0.317290055804, + 0.548033158978, + 0.122794112203 }, + { 1, 1, 1, 1, 2, 2, 3 }, + { 0, 0.75, 1, 2, 0.75, 2, 0.75 }), + power_exp_r( + { 2.16589615432, 1.58417351097, -0.231327054055, 0.0581169164314, + -0.553691372054, 0.489466159094, -0.0242757398435, 0.0624947905017, + -0.121758602252, -0.370556852701, -0.0167758797004, -0.11960736638, + -0.0456193625088, 0.0356127892703, -0.00744277271321, -0.00173957049024, + -0.0218101212895, 0.0243321665592, -0.0374401334235, 0.143387157569, + -0.134919690833, -0.0231512250535, 0.0123631254929, 0.00210583219729, + -0.000339585190264, 0.00559936517716, -0.000303351180556 }, + { 1, 2, 4, 5, 5, 5, 6, 6, 6, 1, 1, 4, 4, 4, 7, 8, 2, 3, 3, 5, 5, 6, 7, 8, 10, 4, 8 }, + { 1.5, 1.5, 2.5, 0, 1.5, 2, 0, 1, 2, 3, 6, 3, 6, 8, + 6, 0, 7, 12, 16, 22, 24, 16, 24, 8, 2, 28, 14 }, + { 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 6 }), + gauss({ -213.654886883, 26641.5691493, -24027.2122046, -283.41603424, 212.472844002 }, + { 2, 2, 2, 3, 3 }, + { 1, 0, 1, 3, 3 }, + { 25, 25, 25, 15, 20 }, + { 1, 1, 1, 1, 1 }, + { 325, 300, 300, 275, 275 }, + { 1.16, 1.19, 1.19, 1.25, 1.22 }), + noan({ -0.666422765408, 0.726086323499, 0.0550686686128 }, + { 3.5, 3.5, 3 }, + { 0.875, 0.925, 0.875 }, + { 0.3, 0.3, 0.3 }, + { 0.7, 0.7, 0.7 }, + { 0.3, 0.3, 1 }, + { 10, 10, 12.5 }, + { 275, 275, 275 }), + + eta_0(0.15178953643112786, + MOLAR_MASS, + 251.196, + 1., + { 0.235156, -0.491266, 0.05211155, 0.05347906, -0.01537102 }), + eta_r({ 1.9036541208525784, + 15.7384720473354, + 0.14207809578440784, + 0.0679058431241662, + -0.030732988514867565 }, + { 0, 0, 3, 0, 1 }, + { 1, 2, 6, 8, 8 }, + { 0, 0, 0, 0, 0 }, + { 1, 1, 1, 1, 0 }), + lambda_r({ 37.0597124660408, + 0.7696647124242399, + 7.5538113451464, + -32.416436589336, + 78.894098855904, + 17.7830586854928, + 107.44756315137599, + 318.39746259479995, + -0.82691726160072, + 2.0846013855224798e-02 }, + { 0.0, 0.0, -1.5, 0.0, -1.0, -1.5, -1.5, -1.5, -3.5, -5.5 }, + { 1.0, 5.0, 1.0, 1.0, 2.0, 0.0, 5.0, 9.0, 0.0, 0.0 }, + { 0, 0, 0, 2, 2, 2, 2, 2, 2, 2 }, + { 0, 0, 0, 5, 5, 5, 5, 5, 5, 5 }) +{ +} + +double +CarbonDioxide::alpha(double delta, double tau) const +{ + // clang-format off + return + this->lead.alpha(delta, tau) + + this->log_tau.alpha(delta, tau) + + this->pe.alpha(delta, tau) + + this->offset.alpha(delta, tau) + + this->power_r.alpha(delta, tau) + + this->power_exp_r.alpha(delta, tau) + + this->gauss.alpha(delta, tau) + + this->noan.alpha(delta, tau); + // clang-format on +} + +double +CarbonDioxide::dalpha_ddelta(double delta, double tau) const +{ + // clang-format off + return + this->lead.ddelta(delta, tau) + + this->power_r.ddelta(delta, tau) + + this->power_exp_r.ddelta(delta, tau) + + this->gauss.ddelta(delta, tau) + + this->noan.ddelta(delta, tau); + // clang-format on +} + +double +CarbonDioxide::dalpha_dtau(double delta, double tau) const +{ + // clang-format off + return + this->lead.dtau(delta, tau) + + this->log_tau.dtau(delta, tau) + + this->pe.dtau(delta, tau) + + this->offset.dtau(delta, tau) + + this->power_r.dtau(delta, tau) + + this->power_exp_r.dtau(delta, tau) + + this->gauss.dtau(delta, tau) + + this->noan.dtau(delta, tau); + // clang-format on +} + +double +CarbonDioxide::d2alpha_ddelta2(double delta, double tau) const +{ + // clang-format off + return + this->lead.d2delta(delta, tau) + + this->power_r.d2delta(delta, tau) + + this->power_exp_r.d2delta(delta, tau) + + this->gauss.d2delta(delta, tau) + + this->noan.d2delta(delta, tau); + // clang-format on +} + +double +CarbonDioxide::d2alpha_dtau2(double delta, double tau) const +{ + // clang-format off + return + this->log_tau.d2tau(delta, tau) + + this->pe.d2tau(delta, tau) + + this->power_r.d2tau(delta, tau) + + this->power_exp_r.d2tau(delta, tau) + + this->gauss.d2tau(delta, tau) + + this->noan.d2tau(delta, tau); + // clang-format on +} + +double +CarbonDioxide::d2alpha_ddeltatau(double delta, double tau) const +{ + // clang-format off + return + this->power_r.d2deltatau(delta, tau) + + this->power_exp_r.d2deltatau(delta, tau) + + this->gauss.d2deltatau(delta, tau) + + this->noan.d2deltatau(delta, tau); + // clang-format on +} + +double +CarbonDioxide::mu_from_rho_T(double rho, double T) const +{ + const double delta = rho / this->rho_c; + const double tau = this->T_c / T; + + double eta = this->eta_0.value(T) + this->eta_r.value(delta, tau); + // [Pa-s] + return eta * 1.0e-6; +} + +double +CarbonDioxide::k_from_rho_T(double rho, double T) const +{ + const double delta = rho / this->rho_c; + const double tau = this->T_c / T; + + double lambda = this->lambda_r.value(delta, tau); + // [W/(m-K)] + return lambda * 1.0e-3; +} + +} // namespace fprops diff --git a/test/src/CMakeLists.txt b/test/src/CMakeLists.txt index 3da9bf6..1d57a1a 100644 --- a/test/src/CMakeLists.txt +++ b/test/src/CMakeLists.txt @@ -1,6 +1,7 @@ add_executable( ${PROJECT_NAME} Air_test.cpp + CarbonDioxide_test.cpp Helium_test.cpp Helmholtz_test.cpp IdealGas_test.cpp diff --git a/test/src/CarbonDioxide_test.cpp b/test/src/CarbonDioxide_test.cpp new file mode 100644 index 0000000..c30878c --- /dev/null +++ b/test/src/CarbonDioxide_test.cpp @@ -0,0 +1,48 @@ +#include "gtest/gtest.h" +#include "CarbonDioxide.h" + +using namespace fprops; + +TEST(CarbonDioxide, p_T) +{ + CarbonDioxide fp; + + double T = 280; + double p = 1e6; + SinglePhaseFluidProperties::Props props = fp.p_T(p, T); + + EXPECT_DOUBLE_EQ(props.p, p); + EXPECT_DOUBLE_EQ(props.T, T); + EXPECT_DOUBLE_EQ(props.rho, 20.199309000812121); + EXPECT_DOUBLE_EQ(props.u, 430888.05580697849); + EXPECT_DOUBLE_EQ(props.cv, 670.91985675762839); + EXPECT_DOUBLE_EQ(props.cp, 925.17988930107481); + EXPECT_DOUBLE_EQ(props.v, 0.049506644012416195); + EXPECT_DOUBLE_EQ(props.s, 2225.7418437948568); + EXPECT_DOUBLE_EQ(props.h, 480394.69981939468); + EXPECT_DOUBLE_EQ(props.w, 252.32654137014907); + EXPECT_DOUBLE_EQ(props.mu, 1.4150532169060509e-05); + EXPECT_DOUBLE_EQ(props.k, 0.015727797767537487); +} + +TEST(CarbonDioxide, v_u) +{ + CarbonDioxide fp; + + double v = 0.049506644012416195; + double u = 430888.05580697849; + SinglePhaseFluidProperties::Props props = fp.v_u(v, u); + + EXPECT_DOUBLE_EQ(props.v, v); + EXPECT_DOUBLE_EQ(props.u, u); + EXPECT_DOUBLE_EQ(props.p, 1e6); + EXPECT_DOUBLE_EQ(props.T, 280); + EXPECT_DOUBLE_EQ(props.rho, 20.199309000812121); + EXPECT_DOUBLE_EQ(props.cv, 670.91985675762839); + EXPECT_DOUBLE_EQ(props.cp, 925.17988930107481); + EXPECT_DOUBLE_EQ(props.s, 2225.7418437948568); + EXPECT_DOUBLE_EQ(props.h, 480394.69981939468); + EXPECT_DOUBLE_EQ(props.w, 252.32654137014907); + EXPECT_DOUBLE_EQ(props.mu, 1.4150532169060509e-05); + EXPECT_DOUBLE_EQ(props.k, 0.015727797767537487); +} From cd14193e42d0ec7444f8af1ad8348f0d5d6b8e5a Mon Sep 17 00:00:00 2001 From: David Andrs Date: Mon, 22 May 2023 15:06:17 -0600 Subject: [PATCH 4/5] Adding carbon dioxide fluid properties (python) --- python/src/fprops.cpp | 6 ++++++ python/src/pyfprops/__init__.py | 1 + python/tests/test_co2.py | 10 ++++++++++ 3 files changed, 17 insertions(+) create mode 100644 python/tests/test_co2.py diff --git a/python/src/fprops.cpp b/python/src/fprops.cpp index 195be19..fa17a3c 100644 --- a/python/src/fprops.cpp +++ b/python/src/fprops.cpp @@ -5,6 +5,7 @@ #include "IdealGas.h" #include "Nitrogen.h" #include "Helium.h" +#include "CarbonDioxide.h" using namespace fprops; @@ -50,4 +51,9 @@ PYBIND11_MODULE(pyfprops, m) .def(py::init()) .def("p_T", &Helium::p_T) .def("v_u", &Helium::v_u); + + py::class_(m, "CarbonDioxide") + .def(py::init()) + .def("p_T", &CarbonDioxide::p_T) + .def("v_u", &CarbonDioxide::v_u); } diff --git a/python/src/pyfprops/__init__.py b/python/src/pyfprops/__init__.py index 7b81c25..8fd732e 100644 --- a/python/src/pyfprops/__init__.py +++ b/python/src/pyfprops/__init__.py @@ -2,6 +2,7 @@ __all__ = [ Air, + CarbonDioxide, Helium, IdealGas, Nitrogen diff --git a/python/tests/test_co2.py b/python/tests/test_co2.py new file mode 100644 index 0000000..76166ee --- /dev/null +++ b/python/tests/test_co2.py @@ -0,0 +1,10 @@ +import pyfprops as fp +import pytest +import math + +def test_co2_valid(): + co2 = fp.CarbonDioxide() + state = co2.p_T(1e6, 280) + assert state.p == 1e6 + assert state.T == 280 + assert math.isclose(state.rho, 20.199309000812121, abs_tol=1e-14) From d55d10a7bcab085462bf4a7733f6db064ec3b95a Mon Sep 17 00:00:00 2001 From: David Andrs Date: Mon, 22 May 2023 15:06:30 -0600 Subject: [PATCH 5/5] Adding carbon dioxide fluid properties (doco) --- docs/CMakeLists.txt | 26 ++++++++++++++++++++- docs/classes/carbon-dioxide.rst | 5 +++++ docs/fluids/carbon-dioxide.rst | 32 ++++++++++++++++++++++++++ docs/pyplots/co2.yml | 17 ++++++++++++++ docs/references.bib | 40 +++++++++++++++++++++++++++++++++ 5 files changed, 119 insertions(+), 1 deletion(-) create mode 100644 docs/classes/carbon-dioxide.rst create mode 100644 docs/fluids/carbon-dioxide.rst create mode 100644 docs/pyplots/co2.yml diff --git a/docs/CMakeLists.txt b/docs/CMakeLists.txt index 481b72f..08e1174 100644 --- a/docs/CMakeLists.txt +++ b/docs/CMakeLists.txt @@ -34,7 +34,6 @@ if(DOXYGEN_FOUND) air_err_u.png air_err_w.png ) - file(GLOB_RECURSE RST_FILES ${PROJECT_SOURCE_DIR}/docs/*.rst) set(HE_ERR_IMGS he_err_c_p.png he_err_c_v.png @@ -46,6 +45,17 @@ if(DOXYGEN_FOUND) he_err_u.png he_err_w.png ) + set(CO2_ERR_IMGS + co2_err_c_p.png + co2_err_c_v.png + co2_err_h.png + co2_err_k.png + co2_err_mu.png + co2_err_rho.png + co2_err_s.png + co2_err_u.png + co2_err_w.png + ) add_custom_command( OUTPUT @@ -59,6 +69,7 @@ if(DOXYGEN_FOUND) ${N2_ERR_IMGS} ${AIR_ERR_IMGS} ${HE_ERR_IMGS} + ${CO2_ERR_IMGS} ) add_custom_command( @@ -110,6 +121,19 @@ if(DOXYGEN_FOUND) ${PROJECT_SOURCE_DIR}/docs/pyplots/he.yml ) + add_custom_command( + OUTPUT + ${CO2_ERR_IMGS} + COMMAND + ${CMAKE_COMMAND} -E env "PYTHONPATH=${PROJECT_BINARY_DIR}/python/src:$ENV{PYTHONPATH}" + python3 plot_err.py co2.yml + WORKING_DIRECTORY + ${PROJECT_SOURCE_DIR}/docs/pyplots + DEPENDS + ${PROJECT_SOURCE_DIR}/docs/pyplots/plot_err.py + ${PROJECT_SOURCE_DIR}/docs/pyplots/co2.yml + ) + add_custom_command( TARGET doc POST_BUILD diff --git a/docs/classes/carbon-dioxide.rst b/docs/classes/carbon-dioxide.rst new file mode 100644 index 0000000..48bc8e9 --- /dev/null +++ b/docs/classes/carbon-dioxide.rst @@ -0,0 +1,5 @@ +Carbon Dioxide +============== + +.. doxygenclass:: fprops::CarbonDioxide + :members: diff --git a/docs/fluids/carbon-dioxide.rst b/docs/fluids/carbon-dioxide.rst new file mode 100644 index 0000000..87d5de3 --- /dev/null +++ b/docs/fluids/carbon-dioxide.rst @@ -0,0 +1,32 @@ +Carbon Dioxide +============== + +Range of validity: + +- Temperature: `216.592 - 1100 K` +- Pressure: to `800 MPa` + +Computation of thermodynamic properties is based on :cite:t:`10.1063/1.555991`. +Viscosity is based on :cite:t:`10.1063/1.556013` and +thermal conductivity is based on :cite:t:`10.1063/1.2213631`. + +Verification against CoolProp package +------------------------------------- + +.. image:: ../pyplots/co2_err_rho.png + +.. image:: ../pyplots/co2_err_u.png + +.. image:: ../pyplots/co2_err_h.png + +.. image:: ../pyplots/co2_err_s.png + +.. image:: ../pyplots/co2_err_w.png + +.. image:: ../pyplots/co2_err_c_p.png + +.. image:: ../pyplots/co2_err_c_v.png + +.. image:: ../pyplots/co2_err_mu.png + +.. image:: ../pyplots/co2_err_k.png diff --git a/docs/pyplots/co2.yml b/docs/pyplots/co2.yml new file mode 100644 index 0000000..f783d26 --- /dev/null +++ b/docs/pyplots/co2.yml @@ -0,0 +1,17 @@ +fprops: + name: co2 + coolprop: CO2 + fprops: CarbonDioxide + p: [ + 0.1, 0.2, + 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, + 5, 6, 8, + 10, 15, 20, 25, + 50, 75, 100, + 200, 400, 600, 800 + ] + T: [ + 220, 230, 240, 250, 260, 270, 280, 290, 300, 310, 320, 330, 340, 350, + 400, 450, 500, 550, + 600, 700, 800, 900, 1000, 1100 + ] diff --git a/docs/references.bib b/docs/references.bib index fc47895..656a3cf 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -57,3 +57,43 @@ @techreport{arp1989thermophysical year={1989}, institution={} } + +@article{10.1063/1.555991, + year = {1996}, + keywords = {CO2}, + title = {{A New Equation of State for Carbon Dioxide Covering the Fluid Region from the Triple‐Point Temperature to 1100 K at Pressures up to 800 MPa}}, + author = {Span, Roland and Wagner, Wolfgang}, + journal = {Journal of Physical and Chemical Reference Data}, + issn = {0047-2689}, + doi = {10.1063/1.555991}, + pages = {1509--1596}, + number = {6}, + volume = {25}, +} + +@article{10.1063/1.556013, + year = {1998}, + keywords = {CO2}, + title = {{The Viscosity of Carbon Dioxide}}, + author = {Fenghour, A and Wakeham, William A and Vesovic, V}, + journal = {Journal of Physical and Chemical Reference Data}, + issn = {0047-2689}, + doi = {10.1063/1.556013}, + pages = {31--44}, + number = {1}, + volume = {27}, +} + +@article{10.1063/1.2213631, + year = {2006}, + keywords = {CO2}, + title = {{A Reference Multiparameter Thermal Conductivity Equation for Carbon Dioxide with an Optimized Functional Form}}, + author = {Scalabrin, G. and Marchi, P. and Finezzo, F. and Span, R.}, + journal = {Journal of Physical and Chemical Reference Data}, + issn = {0047-2689}, + doi = {10.1063/1.2213631}, + abstract = {{A new thermal conductivity equation λ=λ(T,ρ) in a multiparameter format was developed for carbon dioxide through the application of an optimization technique of the functional form. The proposed equation is valid for temperatures from the triple point (Tt=216.592K; Pt=0.51795MPa) to 1000K and pressures up to 200MPa. The calculation of density, which is an independent variable of the equation, from the experimental (T,P) conditions is performed with a high accuracy equation of state for the fluid. The thermal conductivity equation shows an average absolute deviation of 1.19\% on the selected 1407 primary data points. Its performances are slightly better than those of the corresponding conventional model by Vesovic et al. [J. Phys. Chem. Ref. Data 19, 763 (1990)] available from the literature; moreover the new equation is simpler to use in particular for the near-critical region. Tables of generated values of carbon dioxide thermal conductivity are provided for check of the code implementations and for quick evaluations.}}, + pages = {1549--1575}, + number = {4}, + volume = {35}, +}