Add functions for bubble and dew point pressures

prehner · prehner · commit 68db1987bd84 · 2025-01-27T09:28:05.000+01:00
diff --git a/src/lib.rs b/src/lib.rs
@@ -9,4 +9,7 @@ pub use core::{
 #[cfg(feature = "parameter_fit")]
 mod parameter_fit;
 #[cfg(feature = "parameter_fit")]
-pub use parameter_fit::{equilibrium_liquid_density, liquid_density, vapor_pressure};
+pub use parameter_fit::{
+    bubble_point_pressure, dew_point_pressure, equilibrium_liquid_density, liquid_density,
+    vapor_pressure,
+};
diff --git a/src/parameter_fit/mod.rs b/src/parameter_fit/mod.rs
@@ -19,10 +19,10 @@ pub fn vapor_pressure<R: ResidualHelmholtzEnergy<1>, const P: usize>(
     let v2 = 1.0 / vle.vapor().density.to_reduced();
     let t = temperature.into_reduced();
     let (a1, a2) = {
-        let t = DualVec::from(t);
-        let v1 = DualVec::from(v1);
-        let v2 = DualVec::from(v2);
-        let x = SVector::from([DualVec::from(1.0)]);
+        let t = Gradient::from(t);
+        let v1 = Gradient::from(v1);
+        let v2 = Gradient::from(v2);
+        let x = SVector::from([Gradient::from(1.0)]);
 
         let a1 = R::residual_molar_helmholtz_energy(&eos.parameters, t, v1, &x);
         let a2 = R::residual_molar_helmholtz_energy(&eos.parameters, t, v2, &x);
@@ -43,10 +43,10 @@ pub fn equilibrium_liquid_density<R: ResidualHelmholtzEnergy<1>, const P: usize>
     let v_v = 1.0 / vle.vapor().density.to_reduced();
     let t = temperature.into_reduced();
     let (f_l, p_l, dp_l, a_v) = {
-        let t = DualVec::from(temperature.into_reduced());
-        let v_l = DualVec::from(v_l);
-        let v_v = DualVec::from(v_v);
-        let x = SVector::from([DualVec::from(1.0)]);
+        let t = Gradient::from(temperature.into_reduced());
+        let v_l = Gradient::from(v_l);
+        let v_v = Gradient::from(v_v);
+        let x = SVector::from([Gradient::from(1.0)]);
 
         let (f_l, p_l, dp_l) = R::dp_drho(&eos.parameters, t, v_l, &x);
         let a_v = R::residual_molar_helmholtz_energy(&eos.parameters, t, v_v, &x);
@@ -70,9 +70,9 @@ pub fn liquid_density<R: ResidualHelmholtzEnergy<1>, const P: usize>(
     let v = 1.0 / state.density.to_reduced();
     let p0 = pressure.into_reduced();
     let (p, dp) = {
-        let t = DualVec::from(t);
-        let v = DualVec::from(v);
-        let x = SVector::from([DualVec::from(1.0)]);
+        let t = Gradient::from(t);
+        let v = Gradient::from(v);
+        let x = SVector::from([Gradient::from(1.0)]);
         let (_, p, dp) = R::dp_drho(&eos.parameters, t, v, &x);
 
         (p, dp)
@@ -82,11 +82,101 @@ pub fn liquid_density<R: ResidualHelmholtzEnergy<1>, const P: usize>(
     Ok(Density::from_reduced(rho))
 }
 
+pub fn bubble_point_pressure<R: ResidualHelmholtzEnergy<2>, const P: usize>(
+    eos: &HelmholtzEnergyWrapper<R, Gradient<P>, 2>,
+    temperature: Temperature,
+    pressure: Option<Pressure>,
+    liquid_molefracs: SVector<f64, 2>,
+) -> EosResult<Pressure<Gradient<P>>> {
+    let x = arr1(liquid_molefracs.as_slice());
+    let vle = PhaseEquilibrium::bubble_point(
+        &eos.eos,
+        temperature,
+        &x,
+        pressure,
+        None,
+        Default::default(),
+    )?;
+
+    let v_l = 1.0 / vle.liquid().density.to_reduced();
+    let v_v = 1.0 / vle.vapor().density.to_reduced();
+    let y = &vle.vapor().molefracs;
+    let y: SVector<_, 2> = SVector::from_fn(|i, _| y[i]);
+    let t = temperature.into_reduced();
+    let (a_l, a_v, v_l, v_v) = {
+        let t = Gradient::from(t);
+        let v_l = Gradient::from(v_l);
+        let v_v = Gradient::from(v_v);
+        let y = y.map(Gradient::from);
+        let x = liquid_molefracs.map(Gradient::from);
+
+        let a_v = R::residual_molar_helmholtz_energy(&eos.parameters, t, v_v, &y);
+        let (p_l, mu_res_l, dp_l, dmu_l) = R::dmu_dv(&eos.parameters, t, v_l, &x);
+        let vi_l = dmu_l / dp_l;
+        let v_l = vi_l.dot(&y);
+        let a_l = (mu_res_l - vi_l * p_l).dot(&y);
+        (a_l, a_v, v_l, v_v)
+    };
+    let rho_l = vle.liquid().partial_density.to_reduced();
+    let rho_l = [rho_l[0], rho_l[1]];
+    let rho_v = vle.vapor().partial_density.to_reduced();
+    let rho_v = [rho_v[0], rho_v[1]];
+    let p = -(a_v - a_l
+        + t * (y[0] * (rho_v[0] / rho_l[0]).ln() + y[1] * (rho_v[1] / rho_l[1]).ln() - 1.0))
+        / (v_v - v_l);
+    Ok(Pressure::from_reduced(p))
+}
+
+pub fn dew_point_pressure<R: ResidualHelmholtzEnergy<2>, const P: usize>(
+    eos: &HelmholtzEnergyWrapper<R, Gradient<P>, 2>,
+    temperature: Temperature,
+    pressure: Option<Pressure>,
+    vapor_molefracs: SVector<f64, 2>,
+) -> EosResult<Pressure<Gradient<P>>> {
+    let y = arr1(vapor_molefracs.as_slice());
+    let vle = PhaseEquilibrium::dew_point(
+        &eos.eos,
+        temperature,
+        &y,
+        pressure,
+        None,
+        Default::default(),
+    )?;
+
+    let v_l = 1.0 / vle.liquid().density.to_reduced();
+    let v_v = 1.0 / vle.vapor().density.to_reduced();
+    let x = &vle.liquid().molefracs;
+    let x: SVector<_, 2> = SVector::from_fn(|i, _| x[i]);
+    let t = temperature.into_reduced();
+    let (a_l, a_v, v_l, v_v) = {
+        let t = Gradient::from(t);
+        let v_l = Gradient::from(v_l);
+        let v_v = Gradient::from(v_v);
+        let x = x.map(Gradient::from);
+        let y = vapor_molefracs.map(Gradient::from);
+
+        let a_l = R::residual_molar_helmholtz_energy(&eos.parameters, t, v_l, &x);
+        let (p_v, mu_res_v, dp_v, dmu_v) = R::dmu_dv(&eos.parameters, t, v_v, &y);
+        let vi_v = dmu_v / dp_v;
+        let v_v = vi_v.dot(&x);
+        let a_v = (mu_res_v - vi_v * p_v).dot(&x);
+        (a_l, a_v, v_l, v_v)
+    };
+    let rho_l = vle.liquid().partial_density.to_reduced();
+    let rho_l = [rho_l[0], rho_l[1]];
+    let rho_v = vle.vapor().partial_density.to_reduced();
+    let rho_v = [rho_v[0], rho_v[1]];
+    let p = -(a_l - a_v
+        + t * (x[0] * (rho_l[0] / rho_v[0]).ln() + x[1] * (rho_l[1] / rho_v[1]).ln() - 1.0))
+        / (v_l - v_v);
+    Ok(Pressure::from_reduced(p))
+}
+
 #[cfg(test)]
 mod test {
     use super::*;
-    use crate::eos::pcsaft::test::{pcsaft, pcsaft_non_assoc};
-    use crate::eos::PcSaftPure;
+    use crate::eos::pcsaft::test::{pcsaft, pcsaft_binary, pcsaft_non_assoc};
+    use crate::eos::{PcSaftBinary, PcSaftPure};
     use crate::{ParametersAD, PhaseEquilibriumAD, StateAD};
     use approx::assert_relative_eq;
     use nalgebra::U1;
@@ -246,4 +336,64 @@ mod test {
         }
         Ok(())
     }
+
+    #[test]
+    fn test_bubble_point_pressure() -> EosResult<()> {
+        let (pcsaft, _) = pcsaft_binary()?;
+        let pcsaft = pcsaft.wrap();
+        let pcsaft_ad = pcsaft.named_derivatives(["k_ij"]);
+        let temperature = 500.0 * KELVIN;
+        let x = SVector::from([0.5, 0.5]);
+        let p = bubble_point_pressure(&pcsaft_ad, temperature, None, x)?;
+        let p = p.convert_into(BAR);
+        let (p, [[grad]]) = (p.re, p.eps.unwrap_generic(U1, U1).data.0);
+
+        println!("{:.5}", p);
+        println!("{:.5?}", grad);
+
+        let (params, mut kij) = pcsaft.parameters;
+        let h = 1e-7;
+        kij += h;
+        let pcsaft_h = PcSaftBinary::new(params, kij).wrap();
+        let (_, p_h) = PhaseEquilibriumAD::bubble_point(&pcsaft_h, temperature, x)?;
+        let dp_h = (p_h.convert_into(BAR) - p) / h;
+        println!(
+            "k_ij: {:11.5} {:11.5} {:.3e}",
+            dp_h,
+            grad,
+            ((dp_h - grad) / grad).abs()
+        );
+        assert_relative_eq!(grad, dp_h, max_relative = 1e-6);
+        Ok(())
+    }
+
+    #[test]
+    fn test_dew_point_pressure() -> EosResult<()> {
+        let (pcsaft, _) = pcsaft_binary()?;
+        let pcsaft = pcsaft.wrap();
+        let pcsaft_ad = pcsaft.named_derivatives(["k_ij"]);
+        let temperature = 500.0 * KELVIN;
+        let y = SVector::from([0.5, 0.5]);
+        let p = dew_point_pressure(&pcsaft_ad, temperature, None, y)?;
+        let p = p.convert_into(BAR);
+        let (p, [[grad]]) = (p.re, p.eps.unwrap_generic(U1, U1).data.0);
+
+        println!("{:.5}", p);
+        println!("{:.5?}", grad);
+
+        let (params, mut kij) = pcsaft.parameters;
+        let h = 1e-7;
+        kij += h;
+        let pcsaft_h = PcSaftBinary::new(params, kij).wrap();
+        let (_, p_h) = PhaseEquilibriumAD::dew_point(&pcsaft_h, temperature, y)?;
+        let dp_h = (p_h.convert_into(BAR) - p) / h;
+        println!(
+            "k_ij: {:11.5} {:11.5} {:.3e}",
+            dp_h,
+            grad,
+            ((dp_h - grad) / grad).abs()
+        );
+        assert_relative_eq!(grad, dp_h, max_relative = 1e-6);
+        Ok(())
+    }
 }
