Merge pull request #41 from andrsd/rho-p

Adding rho,p API

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$

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;

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$

python/src/fprops.cpp

@@ -36,6 +36,7 @@ PYBIND11_MODULE(pyfprops, m)
     py::class_<Air>(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,25 +45,29 @@ 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);
 
     py::class_<Nitrogen>(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_<Helium>(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_<CarbonDioxide>(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);
 }

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
 {

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
 {

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;

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;

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;

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;

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;

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;