diff --git a/docs/classes/helmholtz.rst b/docs/classes/helmholtz.rst index 93c8e08..bd06a67 100644 --- a/docs/classes/helmholtz.rst +++ b/docs/classes/helmholtz.rst @@ -3,3 +3,30 @@ Helmholtz .. doxygenclass:: fprops::Helmholtz :members: + +Models +------ + +.. doxygenclass:: fprops::Helmholtz::IdealGasLead + :members: + +.. doxygenclass:: fprops::Helmholtz::IdealGasLogTau + :members: + +.. doxygenclass:: fprops::Helmholtz::IdealGasPower + :members: + +.. doxygenclass:: fprops::Helmholtz::IdealPlanckEinsteinGeneralized + :members: + +.. doxygenclass:: fprops::Helmholtz::IdealGasPlanckEinsteinFunctionT + :members: + +.. doxygenclass:: fprops::Helmholtz::ResidualPower + :members: + +.. doxygenclass:: fprops::Helmholtz::ResidualPowerExp + :members: + +.. doxygenclass:: fprops::Helmholtz::ResidualGaussian + :members: diff --git a/docs/classes/transport-models.rst b/docs/classes/transport-models.rst new file mode 100644 index 0000000..05ae8cd --- /dev/null +++ b/docs/classes/transport-models.rst @@ -0,0 +1,11 @@ +TransportModels +=============== + +.. doxygenclass:: fprops::Eta0AndPoly + :members: + +.. doxygenclass:: fprops::LennardJones + :members: + +.. doxygenclass:: fprops::ModifiedBatshinskiHildebrand + :members: diff --git a/docs/conf.py.in b/docs/conf.py.in index 3df1fab..8291815 100644 --- a/docs/conf.py.in +++ b/docs/conf.py.in @@ -3,7 +3,7 @@ import sys import subprocess project = 'fprops' -copyright = '2022, David Andrs' +copyright = '2022-2023, David Andrs' author = 'David Andrs' master_doc = 'index' diff --git a/docs/fluids.rst b/docs/fluids.rst new file mode 100644 index 0000000..8acff60 --- /dev/null +++ b/docs/fluids.rst @@ -0,0 +1,9 @@ +Fluids +====== + +.. toctree:: + :maxdepth: 2 + :hidden: + :glob: + + fluids/* diff --git a/docs/index.rst b/docs/index.rst index 73ab45d..9cd0d59 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -8,17 +8,9 @@ fprops :hidden: install.rst + fluids.rst references.rst -.. toctree:: - :maxdepth: 1 - :titlesonly: - :hidden: - :caption: Fluids - :glob: - - fluids/* - .. toctree:: :maxdepth: 1 :titlesonly: diff --git a/include/Helmholtz.h b/include/Helmholtz.h index 8af4ddc..4f8f595 100644 --- a/include/Helmholtz.h +++ b/include/Helmholtz.h @@ -1,6 +1,10 @@ #pragma once #include "SinglePhaseFluidProperties.h" +#include "Numerics.h" +#include +#include +#include namespace fprops { @@ -98,6 +102,527 @@ class Helmholtz : public SinglePhaseFluidProperties { const double rho_c; /// Critical temperature \f$[K]\f$ const double T_c; + + /// The leading term in the EOS used to set the desired reference state + /// + /// @tparam T The basic data type + /// + /// \f$\alpha = \ln(\delta) + a_1 + a_2 \tau\f$ + template + class IdealGasLead { + public: + IdealGasLead(double a1, double a2) : a1(a1), a2(a2) {} + + T + alpha(T delta, T tau) const + { + return std::log(delta) + this->a1 + this->a2 * tau; + } + + T + ddelta(T delta, T tau) const + { + return 1.0 / delta; + } + + T + dtau(T delta, T tau) const + { + return this->a2; + } + + T + d2delta(T delta, T tau) const + { + return -1.0 / delta / delta; + } + + private: + double a1, a2; + }; + + /// Log(tau) term + /// + /// @tparam T The basic data type + /// + /// \f$\alpha = a_1\ln(\tau)\f$ + template + class IdealGasLogTau { + public: + explicit IdealGasLogTau(double a1) : a1(a1) {} + + T + alpha(T delta, T tau) const + { + return this->a1 * std::log(tau); + } + + T + dtau(T delta, T tau) const + { + return this->a1 / tau; + } + + T + d2tau(T delta, T tau) const + { + return -this->a1 / tau / tau; + } + + private: + double a1; + }; + + /// Sum of powers used for ideal part of \f$\alpha\f$ + /// + /// @tparam T The basic data type + /// + /// \f$\alpha = \displaystyle\sum_{i=0}^n N_i \tau^{t_i}\f$ + template + class IdealGasPower { + public: + IdealGasPower(const std::vector & n, const std::vector & t) : n(n), t(t) + { + assert(n.size() == t.size()); + } + + T + alpha(T delta, T tau) const + { + T sum = 0; + for (std::size_t i = 0; i < this->n.size(); ++i) + sum += this->n[i] * std::pow(tau, this->t[i]); + return sum; + } + + T + dtau(T delta, T tau) const + { + T sum = 0; + for (std::size_t i = 0; i < n.size(); ++i) + sum += this->n[i] * this->t[i] * std::pow(tau, this->t[i] - 1); + return sum; + } + + T + d2tau(T delta, T tau) const + { + T sum = 0; + for (std::size_t i = 0; i < this->n.size(); ++i) + sum += this->n[i] * this->t[i] * (this->t[i] - 1) * std::pow(tau, this->t[i] - 2); + return sum; + } + + private: + std::vector n; + std::vector t; + }; + + /// Generalized Planck-Einstein model + /// + /// @tparam T The basic data type + /// + /// \f$ \alpha = \displaystyle\sum_{i=0}^{n} N_i \log(c_i + d_i \exp(\theta \tau)) \f$ + template + class IdealPlanckEinsteinGeneralized { + public: + IdealPlanckEinsteinGeneralized(const std::vector & n, + const std::vector & theta, + const std::vector & c, + const std::vector & d) : + n(n), + theta(theta), + c(c), + d(d) + { + assert(n.size() == theta.size()); + assert(n.size() == c.size()); + assert(n.size() == d.size()); + } + + T + alpha(T delta, T tau) const + { + T sum = 0; + for (std::size_t i = 0; i < this->n.size(); i++) + sum += + this->n[i] * std::log(this->c[i] + this->d[i] * std::exp(this->theta[i] * tau)); + return sum; + } + + T + dtau(T delta, T tau) const + { + T sum = 0; + for (std::size_t i = 0; i < this->n.size(); ++i) { + double exp_theta_tau = std::exp(this->theta[i] * tau); + sum += this->n[i] * this->theta[i] * this->d[i] * exp_theta_tau / + (this->c[i] + this->d[i] * exp_theta_tau); + } + return sum; + } + + T + d2tau(T delta, T tau) const + { + T sum = 0; + for (std::size_t i = 0; i < this->n.size(); ++i) { + double exp_theta_tau = std::exp(this->theta[i] * tau); + sum += this->n[i] * sqr(this->theta[i]) * this->c[i] * this->d[i] * exp_theta_tau / + sqr(this->c[i] + this->d[i] * exp_theta_tau); + } + return sum; + } + + private: + std::vector n; + std::vector theta; + std::vector c; + std::vector d; + }; + + /// Planck-Einstein model as a function of temperature + /// + /// @tparam T The basic data type + /// + /// \f$ \alpha = \displaystyle\sum_{i=0}^{n} N_i \log(1 - \exp({v_i \over T_{crit}} \tau)) \f$ + template + class IdealGasPlanckEinsteinFunctionT { + public: + IdealGasPlanckEinsteinFunctionT(double T_crit, + const std::vector & n, + const std::vector & v) + { + std::vector theta(n.size(), 0.); + for (std::size_t i = 0; i < v.size(); ++i) + theta[i] = -v[i] / T_crit; + std::vector c(n.size(), 1); + std::vector d(c.size(), -1); + this->generalized = new IdealPlanckEinsteinGeneralized(n, theta, c, d); + } + + T + alpha(T delta, T tau) const + { + return this->generalized->alpha(delta, tau); + } + + T + dtau(T delta, T tau) const + { + return this->generalized->dtau(delta, tau); + } + + T + d2tau(T delta, T tau) const + { + return this->generalized->d2tau(delta, tau); + } + + private: + IdealPlanckEinsteinGeneralized * generalized; + }; + + /// Power model + /// + /// @tparam T basic data type + /// + /// \f$ \alpha = \displaystyle\sum_{i=0}^{n} N_i \delta^{d_i} \tau^{t_i} \f$ + template + class ResidualPower { + public: + ResidualPower(const std::vector & n, + const std::vector & d, + const std::vector & t) : + n(n), + d(d), + t(t) + { + assert(n.size() == d.size()); + assert(n.size() == t.size()); + } + + T + alpha(T delta, T tau) const + { + T sum = 0; + for (std::size_t i = 0; i < this->n.size(); i++) + sum += this->n[i] * std::pow(delta, this->d[i]) * std::pow(tau, this->t[i]); + return sum; + } + + T + ddelta(T delta, T tau) const + { + T sum = 0; + for (std::size_t i = 0; i < this->n.size(); i++) + sum += this->n[i] * this->d[i] * std::pow(delta, this->d[i] - 1) * + std::pow(tau, this->t[i]); + return sum; + } + + T + dtau(T delta, T tau) const + { + T sum = 0; + for (std::size_t i = 0; i < this->n.size(); i++) + sum += this->n[i] * std::pow(delta, this->d[i]) * this->t[i] * + std::pow(tau, this->t[i] - 1); + return sum; + } + + T + d2delta(T delta, T tau) const + { + T sum = 0; + for (std::size_t i = 0; i < this->n.size(); i++) + sum += this->n[i] * this->d[i] * (this->d[i] - 1) * + std::pow(delta, this->d[i] - 2) * std::pow(tau, this->t[i]); + return sum; + } + + T + d2tau(T delta, T tau) const + { + T sum = 0; + for (std::size_t i = 0; i < this->n.size(); i++) + sum += this->n[i] * std::pow(delta, this->d[i]) * this->t[i] * (this->t[i] - 1) * + std::pow(tau, this->t[i] - 2); + return sum; + } + + T + d2deltatau(T delta, T tau) const + { + T sum = 0; + for (std::size_t i = 0; i < this->n.size(); i++) + sum += this->n[i] * this->d[i] * std::pow(delta, this->d[i] - 1) * this->t[i] * + std::pow(tau, this->t[i] - 1); + return sum; + } + + private: + std::vector n; + std::vector d; + std::vector t; + }; + + /// Power model with exponential term + /// + /// @tparam T basic data type + /// @tparam L type of the `l`-coefficient + /// + /// \f$ \alpha = \displaystyle\sum_{i=0}^{n} N_i \delta^{d_i} \tau^{t_i} \exp(-\delta^{l_i}) \f$ + template + class ResidualPowerExp { + public: + ResidualPowerExp(const std::vector & n, + const std::vector & d, + const std::vector & t, + const std::vector & l) : + n(n), + d(d), + t(t), + l(l) + { + assert(n.size() == d.size()); + assert(n.size() == t.size()); + assert(n.size() == l.size()); + } + + T + alpha(T delta, T tau) const + { + T sum = 0; + for (std::size_t i = 0; i < this->n.size(); i++) + sum += this->n[i] * std::pow(delta, this->d[i]) * std::pow(tau, this->t[i]) * + std::exp(-std::pow(delta, this->l[i])); + return sum; + } + + T + ddelta(T delta, T tau) const + { + T sum = 0; + for (std::size_t i = 0; i < this->n.size(); i++) + sum += this->n[i] * std::pow(delta, this->d[i] - 1) * std::pow(tau, this->t[i]) * + std::exp(-std::pow(delta, this->l[i])) * + (this->d[i] - this->l[i] * std::pow(delta, this->l[i])); + return sum; + } + + T + dtau(T delta, T tau) const + { + T sum = 0; + for (std::size_t i = 0; i < this->n.size(); i++) + sum += this->n[i] * this->t[i] * std::pow(delta, this->d[i]) * + std::pow(tau, this->t[i] - 1) * std::exp(-std::pow(delta, this->l[i])); + return sum; + } + + T + d2delta(T delta, T tau) const + { + T sum = 0; + for (std::size_t i = 0; i < this->n.size(); i++) + sum += + this->n[i] * std::pow(delta, this->d[i] - 2) * std::pow(tau, this->t[i]) * + std::exp(-std::pow(delta, this->l[i])) * + (sqr(this->l[i]) * std::pow(delta, 2 * this->l[i]) + + (this->d[i] - 1) * this->d[i] - + this->l[i] * (2 * this->d[i] + this->l[i] - 1) * std::pow(delta, this->l[i])); + return sum; + } + + T + d2tau(T delta, T tau) const + { + T sum = 0; + for (std::size_t i = 0; i < this->n.size(); i++) + sum += this->n[i] * (this->t[i] - 1) * this->t[i] * std::pow(delta, this->d[i]) * + std::exp(-std::pow(delta, this->l[i])) * std::pow(tau, this->t[i] - 2); + return sum; + } + + T + d2deltatau(T delta, T tau) const + { + T sum = 0; + for (std::size_t i = 0; i < this->n.size(); i++) + sum += this->n[i] * this->t[i] * std::pow(delta, this->d[i] - 1) * + std::pow(tau, this->t[i] - 1) * std::exp(-std::pow(delta, this->l[i])) * + (this->d[i] - this->l[i] * std::pow(delta, this->l[i])); + return sum; + } + + private: + std::vector n; + std::vector d; + std::vector t; + std::vector l; + }; + + /// Gaussian model for residual part + /// + /// @tparam T The basic data type + /// + /// \f$ \alpha = \displaystyle\sum_{i=0}^n N_i \delta^{d_i} \tau^{t_i} + /// \exp(-\eta_i (\delta - \epsilon_i)^2 - \beta_i (\tau - \gamma_i)^2) \f$ + template + class ResidualGaussian { + public: + ResidualGaussian(const std::vector & n, + const std::vector & d, + const std::vector & t, + const std::vector & eta, + const std::vector & epsilon, + const std::vector & beta, + const std::vector & gamma) : + n(n), + d(d), + t(t), + eta(eta), + epsilon(epsilon), + beta(beta), + gamma(gamma) + { + assert(n.size() == d.size()); + assert(n.size() == t.size()); + assert(n.size() == eta.size()); + assert(n.size() == epsilon.size()); + assert(n.size() == beta.size()); + assert(n.size() == gamma.size()); + } + + T + alpha(T delta, T tau) const + { + T sum = 0; + for (std::size_t i = 0; i < this->n.size(); i++) + sum += this->n[i] * std::pow(delta, this->d[i]) * std::pow(tau, this->t[i]) * + std::exp(-this->eta[i] * sqr(delta - this->epsilon[i]) - + this->beta[i] * sqr(tau - this->gamma[i])); + return sum; + } + + T + ddelta(T delta, T tau) const + { + T sum = 0; + for (std::size_t i = 0; i < this->n.size(); i++) + sum += this->n[i] * std::pow(delta, this->d[i] - 1) * std::pow(tau, this->t[i]) * + (this->d[i] - 2.0 * delta * this->eta[i] * (delta - this->epsilon[i])) * + std::exp(-this->eta[i] * sqr(delta - this->epsilon[i]) - + this->beta[i] * sqr(tau - this->gamma[i])); + return sum; + } + + T + dtau(T delta, T tau) const + { + T sum = 0; + for (std::size_t i = 0; i < this->n.size(); i++) + sum += this->n[i] * std::pow(delta, this->d[i]) * std::pow(tau, this->t[i] - 1) * + (2.0 * this->beta[i] * (this->gamma[i] - tau) * tau + this->t[i]) * + std::exp(-this->eta[i] * sqr(delta - this->epsilon[i]) - + this->beta[i] * sqr(tau - this->gamma[i])); + return sum; + } + + T + d2delta(T delta, T tau) const + { + T sum = 0; + for (std::size_t i = 0; i < this->n.size(); i++) + sum += this->n[i] * std::pow(delta, this->d[i] - 2) * std::pow(tau, this->t[i]) * + (2 * sqr(delta) * this->eta[i] * + (2 * this->eta[i] * sqr(delta - this->epsilon[i]) - 1) + + sqr(this->d[i]) + + this->d[i] * (4 * delta * this->eta[i] * (this->epsilon[i] - delta) - 1)) * + std::exp(-this->eta[i] * sqr(delta - this->epsilon[i]) - + this->beta[i] * sqr(tau - this->gamma[i])); + return sum; + } + + T + d2tau(T delta, T tau) const + { + T sum = 0; + for (std::size_t i = 0; i < this->n.size(); i++) + sum += this->n[i] * std::pow(delta, this->d[i]) * std::pow(tau, this->t[i] - 2) * + (2 * this->beta[i] * sqr(tau) * + (2 * this->beta[i] * sqr(this->gamma[i] - tau) - 1) + + sqr(this->t[i]) + + this->t[i] * (4 * this->beta[i] * (this->gamma[i] - tau) * tau - 1)) * + std::exp(-this->eta[i] * sqr(delta - this->epsilon[i]) - + this->beta[i] * sqr(tau - this->gamma[i])); + return sum; + } + + T + d2deltatau(T delta, T tau) const + { + T sum = 0; + for (std::size_t i = 0; i < this->n.size(); i++) + sum += this->n[i] * std::pow(delta, this->d[i] - 1) * + std::pow(tau, this->t[i] - 1) * + (2. * this->beta[i] * (this->gamma[i] - tau) * tau + this->t[i]) * + (2. * delta * this->eta[i] * (this->epsilon[i] - delta) + this->d[i]) * + std::exp(-this->beta[i] * sqr(this->gamma[i] - tau) - + this->eta[i] * sqr(delta - this->epsilon[i])); + return sum; + } + + private: + std::vector n; + std::vector d; + std::vector t; + std::vector eta; + std::vector epsilon; + std::vector beta; + std::vector gamma; + }; }; } // namespace fprops diff --git a/include/Nitrogen.h b/include/Nitrogen.h index 0115286..99195ef 100644 --- a/include/Nitrogen.h +++ b/include/Nitrogen.h @@ -1,6 +1,7 @@ #pragma once #include "Helmholtz.h" +#include "TransportModels.h" namespace fprops { @@ -26,6 +27,20 @@ class Nitrogen : public Helmholtz { [[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; + IdealGasPower power_0; + IdealGasPlanckEinsteinFunctionT pefnt; + ResidualPower power_r; + ResidualPowerExp power_exp_r; + ResidualGaussian gauss; + + LennardJones eta_0; + ModifiedBatshinskiHildebrand eta_r; + Eta0AndPoly lambda_0; + ModifiedBatshinskiHildebrand lambda_r; }; } // namespace fprops diff --git a/include/TransportModels.h b/include/TransportModels.h new file mode 100644 index 0000000..d0ea913 --- /dev/null +++ b/include/TransportModels.h @@ -0,0 +1,152 @@ +#pragma once + +#include +#include +#include + +namespace fprops { + +/// Model for computing viscosity (eta_0) +/// +/// @tparam T The basic data type +/// +/// \f$ lambda_0 = A_0 * \eta_0 + \displaystyle\sum_{i=1}^{n} A_i \cdot \tau^{t_i} \f$ +template +class Eta0AndPoly { +public: + /// Eta0 and polynomial model + /// + /// @param A \f$A_i\f$ coefficients + /// @param t \f$t_i\f$ coefficients + /// + /// Note that coefficient \f$t_0\f$ is not used + Eta0AndPoly(const std::vector & A, const std::vector & t) : A(A), t(t) + { + assert(A.size() >= 1); + assert(A.size() == t.size()); + } + + /// Evaluate the model + /// + /// @param eta0 \f$\eta_0\f$ + /// @param tau \f$\tau\f$ + /// @return The computed value + T + value(double eta0, double tau) const + { + T sum = this->A[0] * eta0; + for (unsigned int i = 1; i < A.size(); i++) + sum += this->A[i] * std::pow(tau, this->t[i]); + return sum; + } + +private: + std::vector A; + std::vector t; +}; + +/// Lennard-Jones model for computing viscosity +/// +/// @tparam T The basic data type +/// +/// \f$ \eta_0(T) = \frac{0.0266958 \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 +/// +/// \f$ \Omega(T^*) = \exp(\displaystyle\sum_{i=0}^n b_i \ln(T^*))^i \f$ +/// +/// where \f$T^* = T / (\epsilon / k)) \f$ and \f$\epsilon / k\f$ is the Lennard-Jones +/// energy parameter. +template +class LennardJones { +public: + /// Lennard-Jones model + /// + /// @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) : + M(M), + epsilon_over_k(epsilon_over_k), + sigma(sigma), + b(b) + { + } + + /// Evaluate the model + /// + /// @param temperature Temperature \f$[K]\f$ + /// @return The computed value + T + value(double temperature) const + { + double log_T_star = std::log(temperature / this->epsilon_over_k); + double Omega_T_star = 0; + 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) / + (this->sigma * this->sigma * Omega_T_star); + } + +private: + double M; + double epsilon_over_k; + const double sigma; + std::vector b; +}; + +/// Modified Batshinski-Hildebrand model +/// +/// @tparam T The basic data type +/// +/// \f$ v = \displaystyle\sum_{i=0}^{n} N_i \tau^{t_i} \delta^{d_i} \exp(-\gamma_i \delta^{l_i})\f$ +template +class ModifiedBatshinskiHildebrand { +public: + /// Modified Batshinki-Hildebrand + /// + /// @param n \f$N_i\f$ + /// @param t \f$t_i\f$ + /// @param d \f$d_i\f$ + /// @param gamma \f$\gamma_i\f$ + /// @param l \f$l_i\f$ + ModifiedBatshinskiHildebrand(const std::vector & n, + const std::vector & t, + const std::vector & d, + const std::vector & gamma, + const std::vector & l) : + n(n), + t(t), + d(d), + gamma(gamma), + l(l) + { + } + + /// Evaluate the model + /// + /// @param delta \f$\delta\f$ + /// @param tau \f$\tau\f$ + /// @return Computed value + T + value(double delta, double tau) const + { + double sum = 0.0; + for (unsigned int i = 0; i < n.size(); i++) + sum += this->n[i] * std::pow(tau, this->t[i]) * std::pow(delta, this->d[i]) * + std::exp(-this->gamma[i] * std::pow(delta, this->l[i])); + return sum; + } + +private: + std::vector n; + std::vector t; + std::vector d; + std::vector gamma; + std::vector l; +}; + +} // namespace fprops diff --git a/src/Nitrogen.cpp b/src/Nitrogen.cpp index 7db2259..2d535e2 100644 --- a/src/Nitrogen.cpp +++ b/src/Nitrogen.cpp @@ -1,199 +1,140 @@ #include "Nitrogen.h" -#include "Numerics.h" -#include namespace fprops { -#define len(a) (sizeof(a) / sizeof(a[0])) - -// Coefficients for ideal gas component of the Helmholtz free energy -static const double a[] = { 2.5, -12.76952708, -0.00784163, -1.934819e-4, - -1.247742e-5, 6.678326e-8, 1.012941, 26.65788 }; - -// Coefficients for residual component of the Helmholtz free energy -static const double N1[] = { 0.924803575275, -0.492448489428, 0.661883336938, - -0.192902649201e1, -0.622469309629e-1, 0.349943957581 }; -static const unsigned int i1[] = { 1, 1, 2, 2, 3, 3 }; -static const double j1[] = { 0.25, 0.875, 0.5, 0.875, 0.375, 0.75 }; - -static const double N2[] = { 0.564857472498, -0.161720005987e1, -0.481395031883, - 0.421150636384, -0.161962230825e-1, 0.172100994165, - 0.735448924933e-2, 0.168077305479e-1, -0.107626664179e-2, - -0.137318088513e-1, 0.635466899859e-3, 0.304432279419e-2, - -0.435762336045e-1, -0.723174889316e-1, 0.389644315272e-1, - -0.21220136391e-1, 0.4808822981509e-2, -0.551990017984e-4, - -0.462016716479e-1, -0.300311716011e-2, 0.368825891208e-1, - -0.25585684622e-2, 0.896915264558e-2, -0.44151337035e-2, - 0.133722924858e-2, 0.264832491957e-3 }; -static const unsigned int i2[] = { 1, 1, 1, 3, 3, 4, 6, 6, 7, 7, 8, 8, 1, - 2, 3, 4, 5, 8, 4, 5, 5, 8, 3, 5, 6, 9 }; -static const double j2[] = { 0.5, 0.75, 2.0, 1.25, 3.5, 1.0, 0.5, 3.0, 0.0, - 2.75, 0.75, 2.5, 4.0, 6.0, 6.0, 3.0, 3.0, 6.0, - 16.0, 11.0, 15.0, 12.0, 12.0, 7.0, 4.0, 16.0 }; -static const unsigned int l2[] = { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, - 2, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4 }; - -static const double N3[] = { 0.196688194015e2, - -0.20911560073e2, - 0.167788306989e-1, - 0.262767566274e4 }; -static const unsigned int i3[] = { 1, 1, 3, 2 }; -static const unsigned int j3[] = { 0, 1, 2, 3 }; -static const unsigned int l3[] = { 2, 2, 2, 2 }; -static const double phi3[] = { 20.0, 20.0, 15.0, 25.0 }; -static const double beta3[] = { 325.0, 325.0, 300.0, 275.0 }; -static const double gamma3[] = { 1.16, 1.16, 1.13, 1.25 }; - -/// Coefficients for viscosity -static const double b_mu[] = { 0.431, -0.4623, 0.08406, 0.005341, -0.00331 }; -static const double N_mu[] = { 10.72, 0.03989, 0.001208, -7.402, 4.62 }; -static const double t_mu[] = { 0.1, 0.25, 3.2, 0.9, 0.3 }; -static const double d_mu[] = { 2, 10, 12, 2, 1 }; -static const double l_mu[] = { 0, 1, 1, 2, 3 }; -static const double gamma_mu[] = { 0.0, 1.0, 1.0, 1.0, 1.0 }; - -/// Coefficients for thermal conductivity -static const double N_k[] = { 8.862, 31.11, -73.13, 20.03, -0.7096, 0.2672 }; -static const double t_k[] = { 0.0, 0.03, 0.2, 0.8, 0.6, 1.9 }; -static const unsigned int d_k[] = { 1, 2, 3, 4, 8, 10 }; -static const unsigned int l_k[] = { 0, 0, 1, 2, 2, 2 }; -static const double gamma_k[] = { 0.0, 0.0, 1.0, 1.0, 1.0, 1.0 }; - -Nitrogen::Nitrogen() : Helmholtz(8.314510, 28.01348e-3, 313.299958972, 126.192) {} +static const double MOLAR_MASS = 28.01348e-3; +static const double T_CRIT = 126.192; + +Nitrogen::Nitrogen() : + Helmholtz(8.314510, MOLAR_MASS, 313.299958972, T_CRIT), + lead(-12.76952708, -0.00784163), + log_tau(2.5), + power_0({ -0.0001934819, -1.247742e-05, 6.678326e-08 }, { -1, -2, -3 }), + pefnt(T_CRIT, { 1.012941 }, { 3364.011 }), + power_r({ 0.924803575275, + -0.492448489428, + 0.661883336938, + -1.92902649201, + -0.0622469309629, + 0.349943957581 }, + { 1, 1, 2, 2, 3, 3 }, + { 0.25, 0.875, 0.5, 0.875, 0.375, 0.75 }), + power_exp_r({ 0.564857472498, -1.61720005987, -0.481395031883, 0.421150636384, + -0.0161962230825, 0.172100994165, 0.00735448924933, 0.0168077305479, + -0.00107626664179, -0.0137318088513, 0.000635466899859, 0.00304432279419, + -0.0435762336045, -0.0723174889316, 0.0389644315272, -0.021220136391, + 0.00408822981509, -5.51990017984e-05, -0.0462016716479, -0.00300311716011, + 0.0368825891208, -0.0025585684622, 0.00896915264558, -0.0044151337035, + 0.00133722924858, 0.000264832491957 }, + { 1, 1, 1, 3, 3, 4, 6, 6, 7, 7, 8, 8, 1, 2, 3, 4, 5, 8, 4, 5, 5, 8, 3, 5, 6, 9 }, + { 0.5, 0.75, 2, 1.25, 3.5, 1, 0.5, 3, 0, 2.75, 0.75, 2.5, 4, + 6, 6, 3, 3, 6, 16, 11, 15, 12, 12, 7, 4, 16 }, + { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4 }), + gauss({ 19.6688194015, -20.911560073, 0.0167788306989, 2627.67566274 }, + { 1, 1, 3, 2 }, + { 0, 1, 2, 3 }, + { 20, 20, 15, 25 }, + { 1, 1, 1, 1 }, + { 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_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 }, + { 0.0, 1.0, 1.0, 1.0, 1.0 }, + { 0, 1, 1, 2, 3 }), + lambda_0({ 1.511, 2.117, -3.332 }, { 0, -1, -0.7 }), + lambda_r({ 8.862, 31.11, -73.13, 20.03, -0.7096, 0.2672 }, + { 0.0, 0.03, 0.2, 0.8, 0.6, 1.9 }, + { 1, 2, 3, 4, 8, 10 }, + { 0.0, 0.0, 1.0, 1.0, 1.0, 1.0 }, + { 0, 0, 1, 2, 2, 2 }) +{ +} double Nitrogen::alpha(double delta, double tau) const { - // Ideal gas component of the Helmholtz free energy - double alpha0 = std::log(delta) + a[0] * std::log(tau) + a[1] + a[2] * tau + a[3] / tau + - a[4] / sqr(tau) + a[5] / cb(tau) + a[6] * std::log(1.0 - std::exp(-a[7] * tau)); - - // Residual component of the Helmholtz free energy - double alphar = 0.0; - for (unsigned int i = 0; i < len(N1); i++) - alphar += N1[i] * std::pow(delta, i1[i]) * std::pow(tau, j1[i]); - for (unsigned int i = 0; i < len(N2); i++) - alphar += N2[i] * std::pow(delta, i2[i]) * std::pow(tau, j2[i]) * - std::exp(-std::pow(delta, l2[i])); - for (unsigned int i = 0; i < len(N3); i++) - alphar += N3[i] * std::pow(delta, i3[i]) * std::pow(tau, j3[i]) * - std::exp(-phi3[i] * sqr(delta - 1.0) - beta3[i] * sqr(tau - gamma3[i])); - - return alpha0 + alphar; + // clang-format off + return + this->lead.alpha(delta, tau) + + this->log_tau.alpha(delta, tau) + + this->power_0.alpha(delta, tau) + + this->pefnt.alpha(delta, tau) + + + this->power_r.alpha(delta, tau) + + this->power_exp_r.alpha(delta, tau) + + this->gauss.alpha(delta, tau); + // clang-format on } double Nitrogen::dalpha_ddelta(double delta, double tau) const { - // Ideal gas component of the Helmholtz free energy - double dalpha0 = 1.0 / delta; - - // Residual component of the Helmholtz free energy - double dalphar = 0.0; - for (unsigned int i = 0; i < len(N1); i++) - dalphar += N1[i] * i1[i] * std::pow(delta, i1[i]) * std::pow(tau, j1[i]); - for (unsigned int i = 0; i < len(N2); i++) - dalphar += N2[i] * std::pow(delta, i2[i]) * std::pow(tau, j2[i]) * - std::exp(-std::pow(delta, l2[i])) * (i2[i] - l2[i] * std::pow(delta, l2[i])); - for (unsigned int i = 0; i < len(N3); i++) - dalphar += N3[i] * std::pow(delta, i3[i]) * std::pow(tau, j3[i]) * - std::exp(-phi3[i] * sqr(delta - 1.0) - beta3[i] * sqr(tau - gamma3[i])) * - (i3[i] - 2.0 * delta * phi3[i] * (delta - 1.0)); - - return dalpha0 + dalphar / delta; + // 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); + // clang-format on } double Nitrogen::dalpha_dtau(double delta, double tau) const { - // Ideal gas component of the Helmholtz free energy - const double dalpha0 = a[0] + a[2] * tau - a[3] / tau - 2.0 * a[4] / sqr(tau) - - 3.0 * a[5] / cb(tau) + a[6] * a[7] * tau / (std::exp(a[7] * tau) - 1.0); - - // Residual component of the Helmholtz free energy - double dalphar = 0.0; - for (unsigned int i = 0; i < len(N1); i++) - dalphar += N1[i] * j1[i] * std::pow(delta, i1[i]) * std::pow(tau, j1[i]); - for (unsigned int i = 0; i < len(N2); i++) - dalphar += N2[i] * j2[i] * std::pow(delta, i2[i]) * std::pow(tau, j2[i]) * - std::exp(-std::pow(delta, l2[i])); - for (unsigned int i = 0; i < len(N3); i++) - dalphar += N3[i] * std::pow(delta, i3[i]) * std::pow(tau, j3[i]) * - std::exp(-phi3[i] * sqr(delta - 1.0) - beta3[i] * sqr(tau - gamma3[i])) * - (j3[i] - 2.0 * tau * beta3[i] * (tau - gamma3[i])); - - return (dalpha0 + dalphar) / tau; + // clang-format off + return + this->lead.dtau(delta, tau) + + this->log_tau.dtau(delta, tau) + + this->power_0.dtau(delta, tau) + + this->pefnt.dtau(delta, tau) + + + this->power_r.dtau(delta, tau) + + this->power_exp_r.dtau(delta, tau) + + this->gauss.dtau(delta, tau); + // clang-format on } double Nitrogen::d2alpha_ddelta2(double delta, double tau) const { - // Ideal gas component of the Helmholtz free energy - const double dalpha0 = -1.0 / delta / delta; - - // Residual component of the Helmholtz free energy - double dalphar = 0.0; - for (unsigned int i = 0; i < len(N1); i++) - dalphar += N1[i] * i1[i] * (i1[i] - 1.0) * std::pow(delta, i1[i]) * std::pow(tau, j1[i]); - - for (unsigned int i = 0; i < len(N2); i++) - dalphar += N2[i] * std::pow(delta, i2[i]) * std::pow(tau, j2[i]) * - std::exp(-std::pow(delta, l2[i])) * - ((i2[i] - l2[i] * std::pow(delta, l2[i])) * - (i2[i] - 1.0 - l2[i] * std::pow(delta, l2[i])) - - l2[i] * l2[i] * std::pow(delta, l2[i])); - - for (unsigned int i = 0; i < len(N3); i++) - dalphar += N3[i] * std::pow(delta, i3[i]) * std::pow(tau, j3[i]) * - std::exp(-phi3[i] * sqr(delta - 1.0) - beta3[i] * sqr(tau - gamma3[i])) * - (sqr(i3[i] - 2.0 * delta * phi3[i] * (delta - 1.0)) - i3[i] - - 2.0 * delta * delta * phi3[i]); - - return dalpha0 + dalphar / delta / delta; + // 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); + // clang-format on } double Nitrogen::d2alpha_dtau2(double delta, double tau) const { - // Ideal gas component of the Helmholtz free energy - const double dalpha0 = - -a[0] + 2.0 * a[3] / tau + 6.0 * a[4] / sqr(tau) + 12.0 * a[5] / cb(tau) - - a[6] * a[7] * a[7] * tau * tau * std::exp(a[7] * tau) / sqr(std::exp(a[7] * tau) - 1.0); - - // Residual component of the Helmholtz free energy - double dalphar = 0.0; - for (unsigned int i = 0; i < len(N1); i++) - dalphar += N1[i] * j1[i] * (j1[i] - 1.0) * std::pow(delta, i1[i]) * std::pow(tau, j1[i]); - for (unsigned int i = 0; i < len(N2); i++) - dalphar += N2[i] * j2[i] * (j2[i] - 1.0) * std::pow(delta, i2[i]) * std::pow(tau, j2[i]) * - std::exp(-std::pow(delta, l2[i])); - for (unsigned int i = 0; i < len(N3); i++) - dalphar += N3[i] * std::pow(delta, i3[i]) * std::pow(tau, j3[i]) * - std::exp(-phi3[i] * sqr(delta - 1.0) - beta3[i] * sqr(tau - gamma3[i])) * - (sqr(j3[i] - 2.0 * tau * beta3[i] * (tau - gamma3[i])) - j3[i] - - 2.0 * tau * tau * beta3[i]); - - return (dalpha0 + dalphar) / tau / tau; + // clang-format off + return + this->log_tau.d2tau(delta, tau) + + this->power_0.d2tau(delta, tau) + + this->pefnt.d2tau(delta, tau) + + + this->power_r.d2tau(delta, tau) + + this->power_exp_r.d2tau(delta, tau) + + this->gauss.d2tau(delta, tau); + // clang-format on } double Nitrogen::d2alpha_ddeltatau(double delta, double tau) const { - // Residual component of the Helmholtz free energy (second derivative of ideal - // component wrt delta and tau is 0) - double dalphar = 0.0; - for (unsigned int i = 0; i < len(N1); i++) - dalphar += N1[i] * i1[i] * j1[i] * std::pow(delta, i1[i]) * std::pow(tau, j1[i]); - for (unsigned int i = 0; i < len(N2); i++) - dalphar += N2[i] * j2[i] * std::pow(delta, i2[i]) * std::pow(tau, j2[i]) * - std::exp(-std::pow(delta, l2[i])) * (i2[i] - l2[i] * std::pow(delta, l2[i])); - for (unsigned int i = 0; i < len(N3); i++) - dalphar += N3[i] * std::pow(delta, i3[i]) * std::pow(tau, j3[i]) * - std::exp(-phi3[i] * sqr(delta - 1.0) - beta3[i] * sqr(tau - gamma3[i])) * - (i3[i] - 2.0 * delta * phi3[i] * (delta - 1.0)) * - (j3[i] - 2.0 * tau * beta3[i] * (tau - gamma3[i])); - - return dalphar / delta / tau; + // clang-format off + return + this->power_r.d2deltatau(delta, tau) + + this->power_exp_r.d2deltatau(delta, tau) + + this->gauss.d2deltatau(delta, tau); + // clang-format on } double @@ -201,23 +142,10 @@ Nitrogen::mu_from_rho_T(double rho, double T) const { const double delta = rho / this->rho_c; const double tau = this->T_c / T; - const double logTstar = std::log(T / 98.94); - - // The dilute gas component - double logOmega = 0.0; - for (unsigned int i = 0; i < len(b_mu); i++) - logOmega += b_mu[i] * std::pow(logTstar, i); - const double mu0 = - 0.0266958 * std::sqrt(1000.0 * this->M * T) / (0.3656 * 0.3656 * std::exp(logOmega)); - - // The residual component - double mur = 0.0; - for (unsigned int i = 0; i < len(N_mu); i++) - mur += N_mu[i] * std::pow(tau, t_mu[i]) * std::pow(delta, d_mu[i]) * - std::exp(-gamma_mu[i] * std::pow(delta, l_mu[i])); + double eta = this->eta_0.value(T) + this->eta_r.value(delta, tau); // [Pa-s] - return (mu0 + mur) * 1.0e-6; + return eta * 1.0e-6; } double @@ -226,18 +154,12 @@ Nitrogen::k_from_rho_T(double rho, double T) const const double delta = rho / this->rho_c; const double tau = this->T_c / T; - // The dilute gas component - const double lambda0 = - 1.511 * mu_from_rho_T(0.0, T) * 1.0e6 + 2.117 / tau - 3.332 * std::pow(tau, -0.7); - - // The residual component - double lambdar = 0.0; - for (unsigned int i = 0; i < len(N_k); ++i) - lambdar += N_k[i] * std::pow(tau, t_k[i]) * std::pow(delta, d_k[i]) * - std::exp(-gamma_k[i] * std::pow(delta, l_k[i])); - - // The thermal conductivity (note: critical enhancement not implemented) - return (lambda0 + lambdar) * 1.0e-3; + double eta0 = this->eta_0.value(T); + double lambda = 0; + lambda += this->lambda_0.value(eta0, tau); + 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 1684032..c154164 100644 --- a/test/src/CMakeLists.txt +++ b/test/src/CMakeLists.txt @@ -5,6 +5,7 @@ add_executable( IdealGas_test.cpp Nitrogen_test.cpp Numerics_test.cpp + TransportModels_test.cpp main.cpp ) target_code_coverage(${PROJECT_NAME}) diff --git a/test/src/Helmholtz_test.cpp b/test/src/Helmholtz_test.cpp index 354a92f..ffa0a26 100644 --- a/test/src/Helmholtz_test.cpp +++ b/test/src/Helmholtz_test.cpp @@ -22,17 +22,167 @@ class MockHelmholtz : public Helmholtz { } // namespace -TEST(Helmholtz, p_T_incorrect) +TEST(HelmholtzTest, p_T_incorrect) { MockHelmholtz fp; EXPECT_THROW(auto p = fp.p_T(1e5, -1), std::domain_error); } -TEST(Helmholtz, v_u_incorrect) +TEST(HelmholtzTest, v_u_incorrect) { MockHelmholtz fp; EXPECT_THROW(auto p = fp.v_u(-1, 1), std::domain_error); EXPECT_THROW(auto p = fp.v_u(1, -1), std::domain_error); } + +TEST(HelmholtzTest, ideal_gas_lead) +{ + class Test : public MockHelmholtz { + public: + Test() : MockHelmholtz(), a(10., 100.) + { + double e = std::exp(1.); + EXPECT_DOUBLE_EQ(a.alpha(e, 1), 111.); + EXPECT_DOUBLE_EQ(a.ddelta(10, 0), 0.1); + EXPECT_DOUBLE_EQ(a.dtau(0, 0), 100.); + EXPECT_DOUBLE_EQ(a.d2delta(2, 0), -0.25); + } + + Helmholtz::IdealGasLead a; + }; + + Test t; +} + +TEST(HelmholtzTest, ideal_gas_log_tau) +{ + class Test : public MockHelmholtz { + public: + Test() : MockHelmholtz(), a(11.) + { + double e = std::exp(1.); + EXPECT_DOUBLE_EQ(a.alpha(1., e), 11.); + EXPECT_DOUBLE_EQ(a.dtau(0, 2), 5.5); + EXPECT_DOUBLE_EQ(a.d2tau(0, 2), -2.75); + } + + Helmholtz::IdealGasLogTau a; + }; + + Test t; +} + +TEST(HelmholtzTest, ideal_gas_power) +{ + class Test : public MockHelmholtz { + public: + Test() : MockHelmholtz(), a({ 3, 4 }, { 2, 3 }) + { + EXPECT_DOUBLE_EQ(a.alpha(0., 2), 44.); + EXPECT_DOUBLE_EQ(a.dtau(0, 2), 60); + EXPECT_DOUBLE_EQ(a.d2tau(0, 2), 54); + } + + Helmholtz::IdealGasPower a; + }; + + Test t; +} + +TEST(HelmholtzTest, ideal_gas_planck_einstein_generalized) +{ + class Test : public MockHelmholtz { + public: + Test() : MockHelmholtz(), a({ 2, 3 }, { 3, 2 }, { 3, 4 }, { 5, 6 }) + { + EXPECT_DOUBLE_EQ(a.alpha(0., 1.), 20.9121719702137190); + EXPECT_DOUBLE_EQ(a.dtau(0., 1.), 11.3294239412548978); + EXPECT_DOUBLE_EQ(a.d2tau(0., 1.), 1.41785827094037168); + } + + Helmholtz::IdealPlanckEinsteinGeneralized a; + }; + + Test t; +} + +TEST(HelmholtzTest, ideal_gas_planck_einstein_function_t) +{ + class Test : public MockHelmholtz { + public: + Test() : MockHelmholtz(), a(2, { 2, 3 }, { 3, 4 }) + { + EXPECT_DOUBLE_EQ(a.alpha(0., 1.), -0.941205291457485162); + EXPECT_DOUBLE_EQ(a.dtau(0., 1.), 1.800756606864598643); + EXPECT_DOUBLE_EQ(a.d2tau(0., 1.), -3.835882116252504991); + } + + Helmholtz::IdealGasPlanckEinsteinFunctionT a; + }; + + Test t; +} + +TEST(HelmholtzTest, residual_power) +{ + class Test : public MockHelmholtz { + public: + Test() : MockHelmholtz(), a({ 2, 3 }, { 4, 5 }, { 6, 7 }) + { + EXPECT_DOUBLE_EQ(a.alpha(2., 3.), 233280); + EXPECT_DOUBLE_EQ(a.ddelta(2., 3.), 571536); + EXPECT_DOUBLE_EQ(a.dtau(2., 3.), 536544); + EXPECT_DOUBLE_EQ(a.d2delta(2., 3.), 1119744); + EXPECT_DOUBLE_EQ(a.d2tau(2., 3.), 1057536); + EXPECT_DOUBLE_EQ(a.d2deltatau(2., 3.), 1318032); + } + + Helmholtz::ResidualPower a; + }; + + Test t; +} + +TEST(HelmholtzTest, residual_power_exp) +{ + class Test : public MockHelmholtz { + public: + Test() : MockHelmholtz(), a({ 2, 3 }, { 3, 4 }, { 4, 5 }, { 2, 3 }) + { + EXPECT_DOUBLE_EQ(a.alpha(2., 3.), 27.649904091654395748); + EXPECT_DOUBLE_EQ(a.ddelta(2., 3.), -98.471030918047725031); + EXPECT_DOUBLE_EQ(a.dtau(2., 3.), 38.170817486157493694); + EXPECT_DOUBLE_EQ(a.d2delta(2., 3.), 423.496477990674375469); + EXPECT_DOUBLE_EQ(a.d2tau(2., 3.), 40.344615314965770409); + EXPECT_DOUBLE_EQ(a.d2deltatau(2., 3.), -144.337494863579960335); + } + + Helmholtz::ResidualPowerExp a; + }; + + Test t; +} + +TEST(HelmholtzTest, residual_gaussian) +{ + class Test : public MockHelmholtz { + public: + Test() : + MockHelmholtz(), + a({ 2, 3 }, { 4, 5 }, { 6, 7 }, { 2, 3 }, { 4, 2 }, { 3, 2 }, { 2, 3 }) + { + EXPECT_DOUBLE_EQ(a.alpha(2., 3.), 209952.389617276034850740); + EXPECT_DOUBLE_EQ(a.ddelta(2., 3.), 524883.896172760348507404); + EXPECT_DOUBLE_EQ(a.dtau(2., 3.), 489886.441530895860597038); + EXPECT_DOUBLE_EQ(a.d2delta(2., 3.), -209914.986358776689179657); + EXPECT_DOUBLE_EQ(a.d2tau(2., 3.), 139971.636427909658606910); + EXPECT_DOUBLE_EQ(a.d2deltatau(2., 3.), 1.224704415308958605e6); + } + + Helmholtz::ResidualGaussian a; + }; + + Test t; +} diff --git a/test/src/Nitrogen_test.cpp b/test/src/Nitrogen_test.cpp index 30cf4d6..a89f657 100644 --- a/test/src/Nitrogen_test.cpp +++ b/test/src/Nitrogen_test.cpp @@ -15,16 +15,16 @@ TEST(NitrogenTest, p_T) EXPECT_DOUBLE_EQ(props.p, p); EXPECT_DOUBLE_EQ(props.T, T); - EXPECT_DOUBLE_EQ(props.rho, 12.074993450711256); - EXPECT_DOUBLE_EQ(props.h, 288083.48983915086); - EXPECT_DOUBLE_EQ(props.u, 205267.70993153556); - EXPECT_DOUBLE_EQ(props.s, 6083.1854964551458); - EXPECT_DOUBLE_EQ(props.cp, 1058.6154679415292); - EXPECT_DOUBLE_EQ(props.cv, 745.56949807601234); - EXPECT_DOUBLE_EQ(props.w, 342.35848440003775); - - EXPECT_DOUBLE_EQ(props.mu, 1.709050710929316e-05); - EXPECT_DOUBLE_EQ(props.k, 0.024857419011055305); + EXPECT_DOUBLE_EQ(props.rho, 12.074993451051515); + EXPECT_DOUBLE_EQ(props.h, 288083.48983922758); + EXPECT_DOUBLE_EQ(props.u, 205267.70993394594); + EXPECT_DOUBLE_EQ(props.s, 6083.1854964583363); + EXPECT_DOUBLE_EQ(props.cp, 1058.6154681901673); + EXPECT_DOUBLE_EQ(props.cv, 745.56949823705611); + EXPECT_DOUBLE_EQ(props.w, 342.35848437431741); + + EXPECT_DOUBLE_EQ(props.mu, 1.7090507109297636e-05); + EXPECT_DOUBLE_EQ(props.k, 0.024857419011067187); } TEST(NitrogenTest, v_u) @@ -40,14 +40,14 @@ TEST(NitrogenTest, v_u) EXPECT_DOUBLE_EQ(props.p, p); EXPECT_DOUBLE_EQ(props.T, T); - EXPECT_DOUBLE_EQ(props.rho, 12.074993450711256); - EXPECT_DOUBLE_EQ(props.h, 288083.48983915086); - EXPECT_DOUBLE_EQ(props.u, 205267.70993153556); - EXPECT_DOUBLE_EQ(props.s, 6083.1854964551458); - EXPECT_DOUBLE_EQ(props.cp, 1058.6154679415292); - EXPECT_DOUBLE_EQ(props.cv, 745.56949807601234); - EXPECT_DOUBLE_EQ(props.w, 342.35848440003775); - - EXPECT_DOUBLE_EQ(props.mu, 1.709050710929316e-05); - EXPECT_DOUBLE_EQ(props.k, 0.024857419011055305); + EXPECT_DOUBLE_EQ(props.rho, 12.074993451051515); + EXPECT_DOUBLE_EQ(props.h, 288083.48983922758); + EXPECT_DOUBLE_EQ(props.u, 205267.70993394594); + EXPECT_DOUBLE_EQ(props.s, 6083.1854964583363); + EXPECT_DOUBLE_EQ(props.cp, 1058.6154681901673); + EXPECT_DOUBLE_EQ(props.cv, 745.56949823705611); + EXPECT_DOUBLE_EQ(props.w, 342.35848437431741); + + EXPECT_DOUBLE_EQ(props.mu, 1.7090507109297636e-05); + EXPECT_DOUBLE_EQ(props.k, 0.024857419011067187); } diff --git a/test/src/TransportModels_test.cpp b/test/src/TransportModels_test.cpp new file mode 100644 index 0000000..222f6cd --- /dev/null +++ b/test/src/TransportModels_test.cpp @@ -0,0 +1,22 @@ +#include "gtest/gtest.h" +#include "TransportModels.h" + +using namespace fprops; + +TEST(TransportTest, eta0_and_poly) +{ + Eta0AndPoly eta0p({ 2, 3 }, { 3, 4 }); + EXPECT_DOUBLE_EQ(eta0p.value(12., 2.), 72.); +} + +TEST(TransportTest, lennard_jones) +{ + LennardJones lj(2.e-3, 3., 2., { 2, 3 }); + EXPECT_DOUBLE_EQ(lj.value(2.), 0.0060967411665096916); +} + +TEST(TransportTest, modified_batshinski_hildebrand) +{ + ModifiedBatshinskiHildebrand mbh({ 2, 3 }, { 3, 4 }, { 4, 5 }, { 4, 3 }, { 2, 3 }); + EXPECT_DOUBLE_EQ(mbh.value(2., 3.), 0.000097523945419603); +}