feos_core/ad/properties/
bubble_point_pressure.rs

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

feos_core/ad/properties/
bubble_point_pressure.rs
