Skip to main content

feos_core/ad/properties/
dew_point_pressure.rs

1use super::PropertyAD;
2use crate::ad::Gradient;
3use crate::{Composition, Contributions, FeosResult, PhaseEquilibrium, ReferenceSystem, Residual};
4use nalgebra::allocator::Allocator;
5use nalgebra::{DefaultAllocator, U1};
6use num_dual::{DualNum, DualStruct, Gradients};
7use quantity::{_Pressure, KELVIN, PASCAL, Pressure, Temperature};
8
9/// Dew point pressure of a binary mixture as function of temperature and
10/// molefracs of the first component.
11///
12/// An initial value for the pressure can be passed as optional argument to
13/// increase robustness and speed.
14pub struct DewPointPressure(pub Temperature, pub f64, pub Option<Pressure>);
15
16impl<'a> From<&'a [f64]> for DewPointPressure {
17    fn from(value: &'a [f64]) -> Self {
18        Self(value[0] * KELVIN, value[1], Some(value[2] * PASCAL))
19    }
20}
21
22impl<N: Gradients> PropertyAD<N> for DewPointPressure
23where
24    DefaultAllocator: Allocator<N> + Allocator<U1, N> + Allocator<N, N>,
25    f64: Composition<f64, N>,
26{
27    type Unit = _Pressure;
28    const REFERENCE: Pressure = PASCAL;
29
30    fn evaluate<E: Residual<N, D>, D: DualNum<f64, Inner = f64> + Copy>(
31        &self,
32        eos: &E,
33    ) -> FeosResult<Pressure<D>>
34    where
35        DefaultAllocator: Allocator<N> + Allocator<U1, N> + Allocator<N, N>,
36    {
37        let t = Temperature::from_inner(&self.0);
38        let p = Option::from_inner(&self.2);
39        let (y, _) = self.1.into_molefracs(&eos.re())?;
40        let y = y.map(D::from);
41        let vle = PhaseEquilibrium::dew_point(eos, t, y, p, None, Default::default())?;
42        Ok(vle.vapor().pressure(Contributions::Total))
43    }
44
45    fn evaluate_gradient<E: Residual<N, Gradient<P>>, const P: usize>(
46        &self,
47        eos: &E,
48    ) -> FeosResult<quantity::Quantity<Gradient<P>, Self::Unit>>
49    where
50        DefaultAllocator: Allocator<N> + Allocator<U1, N> + Allocator<N, N>,
51    {
52        let eos_f64 = eos.re();
53        let (vapor_molefracs, _) = self.1.into_molefracs(&eos_f64)?;
54        let vle = PhaseEquilibrium::dew_point(
55            &eos_f64,
56            self.0,
57            &vapor_molefracs,
58            self.2,
59            None,
60            Default::default(),
61        )?;
62
63        let v_l = 1.0 / vle.liquid().density.to_reduced();
64        let v_v = 1.0 / vle.vapor().density.to_reduced();
65        let x = &vle.liquid().molefracs;
66        let t = self.0.into_reduced();
67        let (a_l, a_v, v_l, v_v) = {
68            let t = Gradient::from(t);
69            let v_l = Gradient::from(v_l);
70            let v_v = Gradient::from(v_v);
71            let x = x.map(Gradient::from);
72            let y = vapor_molefracs.map(Gradient::from);
73
74            let a_l = eos.residual_helmholtz_energy(t, v_l, &x);
75            let (p_v, mu_res_v, dp_v, dmu_v) = eos.dmu_dv(t, v_v, &y);
76            let vi_v = dmu_v / dp_v;
77            let v_v = vi_v.dot(&x);
78            let a_v = (mu_res_v - vi_v * p_v).dot(&x);
79            (a_l, a_v, v_l, v_v)
80        };
81        let rho_l = vle.liquid().partial_density().to_reduced();
82        let rho_l = [rho_l[0], rho_l[1]];
83        let rho_v = vle.vapor().partial_density().to_reduced();
84        let rho_v = [rho_v[0], rho_v[1]];
85        let p = -(a_l - a_v
86            + t * (x[0] * (rho_l[0] / rho_v[0]).ln() + x[1] * (rho_l[1] / rho_v[1]).ln() - 1.0))
87            / (v_l - v_v);
88        Ok(Pressure::from_reduced(p))
89    }
90}