feos_core/ad/properties/
vapor_pressure.rs

1use super::PropertyAD;
2use crate::ad::Gradient;
3use crate::{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/// Vapor pressure of a pure component as function of temperature.
10pub struct VaporPressure(pub Temperature);
11
12impl<'a> From<&'a [f64]> for VaporPressure {
13    fn from(value: &'a [f64]) -> Self {
14        Self(value[0] * KELVIN)
15    }
16}
17
18impl<N: Gradients> PropertyAD<N> for VaporPressure
19where
20    DefaultAllocator: Allocator<N> + Allocator<U1, N> + Allocator<N, N>,
21{
22    type Unit = _Pressure;
23    const REFERENCE: Pressure = PASCAL;
24
25    fn evaluate<E: Residual<N, D>, D: DualNum<f64, Inner = f64> + Copy>(
26        &self,
27        eos: &E,
28    ) -> FeosResult<Pressure<D>> {
29        let t = Temperature::from_inner(&self.0);
30        PhaseEquilibrium::pure_t(eos, t, None, Default::default()).map(|(p, _)| p)
31    }
32
33    fn evaluate_gradient<E: Residual<N, Gradient<P>>, const P: usize>(
34        &self,
35        eos: &E,
36    ) -> FeosResult<Pressure<Gradient<P>>> {
37        let eos_f64 = eos.re();
38        let (_, [vapor_density, liquid_density]) =
39            PhaseEquilibrium::pure_t(&eos_f64, self.0, None, Default::default())?;
40
41        // implicit differentiation is implemented here instead of just calling pure_t with dual
42        // numbers, because for the first derivative, we can avoid calculating density derivatives.
43        let v1 = 1.0 / liquid_density.to_reduced();
44        let v2 = 1.0 / vapor_density.to_reduced();
45        let t = self.0.into_reduced();
46        let (a1, a2) = {
47            let t = Gradient::from(t);
48            let v1 = Gradient::from(v1);
49            let v2 = Gradient::from(v2);
50            let x = E::pure_molefracs();
51
52            let a1 = eos.residual_helmholtz_energy(t, v1, &x);
53            let a2 = eos.residual_helmholtz_energy(t, v2, &x);
54            (a1, a2)
55        };
56
57        let p = -(a1 - a2 + t * (v2 / v1).ln()) / (v1 - v2);
58        Ok(Pressure::from_reduced(p))
59    }
60}
feos_core/ad/properties/vapor_pressure.rs

feos_core/ad/properties/
vapor_pressure.rs
