From ba6f2bb2b26e15827eea015ad5a871cdaf9cfa56 Mon Sep 17 00:00:00 2001 From: David Andrs Date: Wed, 23 Aug 2023 17:40:19 -0600 Subject: [PATCH] Adding rho,p API --- include/Helmholtz.h | 8 ++++ include/IdealGas.h | 1 + include/SinglePhaseFluidProperties.h | 6 +++ python/src/fprops.cpp | 5 ++ src/Helmholtz.cpp | 71 ++++++++++++++++++++++++++++ src/IdealGas.cpp | 25 ++++++++++ test/src/Air_test.cpp | 22 +++++++++ test/src/CarbonDioxide_test.cpp | 22 +++++++++ test/src/Helium_test.cpp | 22 +++++++++ test/src/Helmholtz_test.cpp | 7 +++ test/src/IdealGas_test.cpp | 35 ++++++++++++++ test/src/Nitrogen_test.cpp | 22 +++++++++ 12 files changed, 246 insertions(+) diff --git a/include/Helmholtz.h b/include/Helmholtz.h index 26b9657..fb2d273 100644 --- a/include/Helmholtz.h +++ b/include/Helmholtz.h @@ -21,6 +21,7 @@ class Helmholtz : public SinglePhaseFluidProperties { Helmholtz(double R, double M, double rho_c, double T_c); [[nodiscard]] Props rho_T(double rho, double T) const override; + [[nodiscard]] Props rho_p(double rho, double p) const override; [[nodiscard]] Props p_T(double p, double T) const override; [[nodiscard]] Props v_u(double v, double u) const override; [[nodiscard]] Props h_s(double h, double s) const override; @@ -75,6 +76,13 @@ class Helmholtz : public SinglePhaseFluidProperties { /// @return Density \f$[kg/m^3]\f$ [[nodiscard]] double rho_from_p_T(double p, double T) const; + /// Temperature given density and pressure + /// + /// @param rho Density \f$[kg/m^3]\f$ + /// @param p Pressure \f$[Pa]\f$ + /// @return Temperature \f$[K]\f$ + [[nodiscard]] double T_from_rho_p(double rho, double p) const; + /// Tau given specific volume and internal energy /// /// @param v Specific volume \f$[m^3/kg]\f$ diff --git a/include/IdealGas.h b/include/IdealGas.h index 920b2ed..60192dd 100644 --- a/include/IdealGas.h +++ b/include/IdealGas.h @@ -25,6 +25,7 @@ class IdealGas : public SinglePhaseFluidProperties { void set_k(double k); [[nodiscard]] Props rho_T(double rho, double T) const override; + [[nodiscard]] Props rho_p(double rho, double p) const override; [[nodiscard]] Props p_T(double p, double T) const override; [[nodiscard]] Props v_u(double v, double u) const override; [[nodiscard]] Props h_s(double h, double s) const override; diff --git a/include/SinglePhaseFluidProperties.h b/include/SinglePhaseFluidProperties.h index c05849a..bf9646d 100644 --- a/include/SinglePhaseFluidProperties.h +++ b/include/SinglePhaseFluidProperties.h @@ -73,6 +73,12 @@ class SinglePhaseFluidProperties : public FluidProperties { /// @param T Temperature \f$[K]\f$ [[nodiscard]] virtual Props rho_T(double rho, double T) const = 0; + /// Compute thermodynamical state given density and pressure + /// + /// @param rho Density \f$[kg/m^3]\f$ + /// @param p Pressure \f$[Pa]\f$ + [[nodiscard]] virtual Props rho_p(double rho, double p) const = 0; + /// Compute thermodynamical state given pressure and temperature /// /// @param p Pressure \f$[Pa]\f$ diff --git a/python/src/fprops.cpp b/python/src/fprops.cpp index 06aa03a..d5c0851 100644 --- a/python/src/fprops.cpp +++ b/python/src/fprops.cpp @@ -36,6 +36,7 @@ PYBIND11_MODULE(pyfprops, m) py::class_(m, "Air") .def(py::init()) .def("rho_T", &Air::rho_T) + .def("rho_p", &Air::rho_p) .def("p_T", &Air::p_T) .def("v_u", &Air::v_u); @@ -44,6 +45,7 @@ PYBIND11_MODULE(pyfprops, m) .def("set_mu", &IdealGas::set_mu) .def("set_k", &IdealGas::set_k) .def("rho_T", &IdealGas::rho_T) + .def("rho_p", &IdealGas::rho_p) .def("p_T", &IdealGas::p_T) .def("v_u", &IdealGas::v_u) .def("h_s", &IdealGas::h_s); @@ -51,18 +53,21 @@ PYBIND11_MODULE(pyfprops, m) py::class_(m, "Nitrogen") .def(py::init()) .def("rho_T", &Nitrogen::rho_T) + .def("rho_p", &Nitrogen::rho_p) .def("p_T", &Nitrogen::p_T) .def("v_u", &Nitrogen::v_u); py::class_(m, "Helium") .def(py::init()) .def("rho_T", &Helium::rho_T) + .def("rho_p", &Helium::rho_p) .def("p_T", &Helium::p_T) .def("v_u", &Helium::v_u); py::class_(m, "CarbonDioxide") .def(py::init()) .def("rho_T", &CarbonDioxide::rho_T) + .def("rho_p", &CarbonDioxide::rho_p) .def("p_T", &CarbonDioxide::p_T) .def("v_u", &CarbonDioxide::v_u); } diff --git a/src/Helmholtz.cpp b/src/Helmholtz.cpp index 8cadfb3..b4e33ca 100644 --- a/src/Helmholtz.cpp +++ b/src/Helmholtz.cpp @@ -62,6 +62,53 @@ Helmholtz::rho_T(double rho, double T) const return props; } +SinglePhaseFluidProperties::Props +Helmholtz::rho_p(double rho, double p) const +{ + if (rho < 0) + throw std::domain_error("Negative density"); + + const double T = T_from_rho_p(rho, p); + + const double delta = rho / this->rho_c; + const double tau = this->T_c / T; + + const double a = alpha(delta, tau); + const double da_dd = dalpha_ddelta(delta, tau); + const double da_dt = dalpha_dtau(delta, tau); + const double d2a_dt2 = d2alpha_dtau2(delta, tau); + const double d2a_dd2 = d2alpha_ddelta2(delta, tau); + const double d2a_ddt = d2alpha_ddeltatau(delta, tau); + + Props props; + props.rho = rho; + props.v = 1. / rho; + props.p = p; + props.T = T; + // u + props.u = this->R * T * tau * da_dt / this->M; + // h + props.h = this->R * props.T * (tau * da_dt + delta * da_dd) / this->M; + // w + const double n = 2.0 * delta * da_dd + delta * delta * d2a_dd2 - + sqr(delta * da_dd - delta * tau * d2a_ddt) / (tau * tau * d2a_dt2); + props.w = std::sqrt(this->R * props.T * n / this->M); + // cp = dh/dt + props.cp = this->R * + (-tau * tau * d2a_dt2 + sqr(delta * da_dd - delta * tau * d2a_ddt) / + (2.0 * delta * da_dd + delta * delta * d2a_dd2)) / + this->M; + // cv = du/dt + props.cv = -this->R * tau * tau * d2a_dt2 / this->M; + // s = ... + props.s = this->R * (tau * da_dt - a) / this->M; + // mu + props.mu = mu_from_rho_T(rho, props.T); + // k + props.k = k_from_rho_T(rho, props.T); + return props; +} + SinglePhaseFluidProperties::Props Helmholtz::p_T(double p, double T) const { @@ -191,6 +238,30 @@ Helmholtz::rho_from_p_T(double p, double T) const return newton::root(1.0e-2, f, df); } +double +Helmholtz::T_from_rho_p(double rho, double p) const +{ + auto f = [&rho, &p, this](double T) { + const double delta = rho / this->rho_c; + const double tau = this->T_c / T; + + return this->R * rho * T * delta * dalpha_ddelta(delta, tau) / this->M - p; + }; + auto df = [&rho, this](double T) { + const double delta = rho / this->rho_c; + const double tau = this->T_c / T; + const double dtau_dT = - this->T_c / T / T; + + double K = this->R * rho * delta / this->M; + double t1 = dalpha_ddelta(delta, tau); + double t2 = T * d2alpha_ddeltatau(delta, tau) * dtau_dT; + + return K * (t1 + t2); + }; + + return newton::root(275., f, df); +} + double Helmholtz::tau_from_v_u(double v, double u) const { diff --git a/src/IdealGas.cpp b/src/IdealGas.cpp index ce1d23e..1b38df2 100644 --- a/src/IdealGas.cpp +++ b/src/IdealGas.cpp @@ -46,6 +46,31 @@ IdealGas::rho_T(double rho, double T) const return props; } +SinglePhaseFluidProperties::Props +IdealGas::rho_p(double rho, double p) const +{ + if (rho < 0) + throw std::domain_error("Negative density"); + + Props props; + props.rho = rho; + props.p = p; + props.cp = this->cp; + props.cv = this->cv; + props.mu = this->mu; + props.k = this->k; + props.T = p * this->molar_mass / (rho * R); + props.u = this->cv * props.T; + props.v = 1. / props.rho; + const double n = std::pow(props.T, this->gamma) / std::pow(props.p, this->gamma - 1.0); + if (n <= 0) + throw std::domain_error("Invalid log base for computing entropy"); + props.s = this->cv * std::log(n); + props.h = this->cp * props.T; + props.w = std::sqrt(this->cp * R * props.T / (this->cv * this->molar_mass)); + return props; +} + SinglePhaseFluidProperties::Props IdealGas::p_T(double p, double T) const { diff --git a/test/src/Air_test.cpp b/test/src/Air_test.cpp index 08bd558..1e7e06a 100644 --- a/test/src/Air_test.cpp +++ b/test/src/Air_test.cpp @@ -43,6 +43,28 @@ TEST(Air, rho_T) EXPECT_DOUBLE_EQ(props.w, gold1.w); } +TEST(Air, rho_p) +{ + Air fp; + + double rho = 1.1769510785919943; + double p = 101325; + SinglePhaseFluidProperties::Props props = fp.rho_p(rho, p); + + EXPECT_DOUBLE_EQ(props.rho, gold1.rho); + EXPECT_DOUBLE_EQ(props.T, gold1.T); + EXPECT_DOUBLE_EQ(props.p, gold1.p); + EXPECT_DOUBLE_EQ(props.u, gold1.u); + EXPECT_DOUBLE_EQ(props.cv, gold1.cv); + EXPECT_DOUBLE_EQ(props.cp, gold1.cp); + EXPECT_DOUBLE_EQ(props.mu, gold1.mu); + EXPECT_DOUBLE_EQ(props.k, gold1.k); + EXPECT_DOUBLE_EQ(props.v, gold1.v); + EXPECT_DOUBLE_EQ(props.s, gold1.s); + EXPECT_DOUBLE_EQ(props.h, gold1.h); + EXPECT_DOUBLE_EQ(props.w, gold1.w); +} + TEST(Air, p_T) { Air fp; diff --git a/test/src/CarbonDioxide_test.cpp b/test/src/CarbonDioxide_test.cpp index f14bb6d..1abd2ab 100644 --- a/test/src/CarbonDioxide_test.cpp +++ b/test/src/CarbonDioxide_test.cpp @@ -43,6 +43,28 @@ TEST(CarbonDioxide, rho_T) EXPECT_DOUBLE_EQ(props.w, gold1.w); } +TEST(CarbonDioxide, rho_p) +{ + CarbonDioxide fp; + + double rho = 20.199309000812121; + double p = 1e6; + SinglePhaseFluidProperties::Props props = fp.rho_p(rho, p); + + EXPECT_DOUBLE_EQ(props.rho, gold1.rho); + EXPECT_DOUBLE_EQ(props.T, gold1.T); + EXPECT_DOUBLE_EQ(props.p, gold1.p); + EXPECT_DOUBLE_EQ(props.u, gold1.u); + EXPECT_DOUBLE_EQ(props.cv, gold1.cv); + EXPECT_DOUBLE_EQ(props.cp, gold1.cp); + EXPECT_DOUBLE_EQ(props.mu, gold1.mu); + EXPECT_DOUBLE_EQ(props.k, gold1.k); + EXPECT_DOUBLE_EQ(props.v, gold1.v); + EXPECT_DOUBLE_EQ(props.s, gold1.s); + EXPECT_DOUBLE_EQ(props.h, gold1.h); + EXPECT_DOUBLE_EQ(props.w, gold1.w); +} + TEST(CarbonDioxide, p_T) { CarbonDioxide fp; diff --git a/test/src/Helium_test.cpp b/test/src/Helium_test.cpp index 06ea81c..83c124d 100644 --- a/test/src/Helium_test.cpp +++ b/test/src/Helium_test.cpp @@ -43,6 +43,28 @@ TEST(HeliumTest, rho_T) EXPECT_DOUBLE_EQ(props.w, gold1.w); } +TEST(HeliumTest, rho_p) +{ + Helium fp; + + double rho = 1.7109055009783694; + double p = 1.e6; + SinglePhaseFluidProperties::Props props = fp.rho_p(rho, p); + + EXPECT_DOUBLE_EQ(props.rho, gold1.rho); + EXPECT_DOUBLE_EQ(props.T, gold1.T); + EXPECT_DOUBLE_EQ(props.p, gold1.p); + EXPECT_DOUBLE_EQ(props.u, gold1.u); + EXPECT_DOUBLE_EQ(props.cv, gold1.cv); + EXPECT_DOUBLE_EQ(props.cp, gold1.cp); + EXPECT_DOUBLE_EQ(props.mu, gold1.mu); + EXPECT_DOUBLE_EQ(props.k, gold1.k); + EXPECT_DOUBLE_EQ(props.v, gold1.v); + EXPECT_DOUBLE_EQ(props.s, gold1.s); + EXPECT_DOUBLE_EQ(props.h, gold1.h); + EXPECT_DOUBLE_EQ(props.w, gold1.w); +} + TEST(HeliumTest, p_T) { Helium fp; diff --git a/test/src/Helmholtz_test.cpp b/test/src/Helmholtz_test.cpp index 47e34d7..2add302 100644 --- a/test/src/Helmholtz_test.cpp +++ b/test/src/Helmholtz_test.cpp @@ -30,6 +30,13 @@ TEST(HelmholtzTest, rho_T_incorrect) EXPECT_THROW(auto p = fp.rho_T(1, -1), std::domain_error); } +TEST(HelmholtzTest, rho_p_incorrect) +{ + MockHelmholtz fp; + + EXPECT_THROW(auto p = fp.rho_p(-1, 300), std::domain_error); +} + TEST(HelmholtzTest, h_s) { MockHelmholtz fp; diff --git a/test/src/IdealGas_test.cpp b/test/src/IdealGas_test.cpp index 21b544f..8402a82 100644 --- a/test/src/IdealGas_test.cpp +++ b/test/src/IdealGas_test.cpp @@ -46,6 +46,41 @@ TEST(IdealGas, rho_T) EXPECT_DOUBLE_EQ(props.w, gold1.w); } +TEST(IdealGas, rho_p) +{ + double gamma = 1.4; + double molar_mass = 29.0e-3; + IdealGas fp(gamma, molar_mass); + fp.set_mu(18.23e-6); + fp.set_k(25.68e-3); + + double rho = 0.89892258591830565; + double p = 101325; + SinglePhaseFluidProperties::Props props = fp.rho_p(rho, p); + + EXPECT_DOUBLE_EQ(props.rho, gold1.rho); + EXPECT_DOUBLE_EQ(props.T, gold1.T); + EXPECT_DOUBLE_EQ(props.p, gold1.p); + EXPECT_DOUBLE_EQ(props.u, gold1.u); + EXPECT_DOUBLE_EQ(props.cv, gold1.cv); + EXPECT_DOUBLE_EQ(props.cp, gold1.cp); + EXPECT_DOUBLE_EQ(props.mu, gold1.mu); + EXPECT_DOUBLE_EQ(props.k, gold1.k); + EXPECT_DOUBLE_EQ(props.v, gold1.v); + EXPECT_DOUBLE_EQ(props.s, gold1.s); + EXPECT_DOUBLE_EQ(props.h, gold1.h); + EXPECT_DOUBLE_EQ(props.w, gold1.w); +} + +TEST(IdealGas, rho_p_incorrect) +{ + double gamma = 1.4; + double molar_mass = 29.0e-3; + IdealGas fp(gamma, molar_mass); + + EXPECT_THROW(auto p = fp.rho_p(-1, 1e5), std::domain_error); +} + TEST(IdealGas, p_T) { double gamma = 1.4; diff --git a/test/src/Nitrogen_test.cpp b/test/src/Nitrogen_test.cpp index f47a9fd..e4a5d7a 100644 --- a/test/src/Nitrogen_test.cpp +++ b/test/src/Nitrogen_test.cpp @@ -43,6 +43,28 @@ TEST(NitrogenTest, rho_T) EXPECT_DOUBLE_EQ(props.w, gold1.w); } +TEST(NitrogenTest, rho_p) +{ + Nitrogen fp; + + double rho = 12.074993451051515; + double p = 1.0e6; + SinglePhaseFluidProperties::Props props = fp.rho_p(rho, p); + + EXPECT_DOUBLE_EQ(props.rho, gold1.rho); + EXPECT_DOUBLE_EQ(props.T, gold1.T); + EXPECT_DOUBLE_EQ(props.p, gold1.p); + EXPECT_DOUBLE_EQ(props.u, gold1.u); + EXPECT_DOUBLE_EQ(props.cv, gold1.cv); + EXPECT_DOUBLE_EQ(props.cp, gold1.cp); + EXPECT_DOUBLE_EQ(props.mu, gold1.mu); + EXPECT_DOUBLE_EQ(props.k, gold1.k); + EXPECT_NEAR(props.v, gold1.v, 1e-9); + EXPECT_DOUBLE_EQ(props.s, gold1.s); + EXPECT_DOUBLE_EQ(props.h, gold1.h); + EXPECT_DOUBLE_EQ(props.w, gold1.w); +} + TEST(NitrogenTest, p_T) { Nitrogen fp;