Skip to main content

feos_core/ad/properties/
boiling_temperature.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, first_derivative, partial2};
7use quantity::{_Temperature, KELVIN, PASCAL, Pressure, Temperature};
8
9/// Boiling temperature of a pure component as function of pressure.
10pub struct BoilingTemperature(pub Pressure);
11
12impl<'a> From<&'a [f64]> for BoilingTemperature {
13    fn from(value: &'a [f64]) -> Self {
14        Self(value[0] * PASCAL)
15    }
16}
17
18impl<N: Gradients> PropertyAD<N> for BoilingTemperature
19where
20    DefaultAllocator: Allocator<N> + Allocator<U1, N> + Allocator<N, N>,
21{
22    type Unit = _Temperature;
23    const REFERENCE: Temperature = KELVIN;
24
25    fn evaluate<E: Residual<N, D>, D: DualNum<f64, Inner = f64> + Copy>(
26        &self,
27        eos: &E,
28    ) -> FeosResult<Temperature<D>> {
29        let p = Pressure::from_inner(&self.0);
30        PhaseEquilibrium::pure_p(eos, p, None, Default::default()).map(|(t, _)| t)
31    }
32
33    fn evaluate_gradient<E: Residual<N, Gradient<P>>, const P: usize>(
34        &self,
35        eos: &E,
36    ) -> FeosResult<quantity::Quantity<Gradient<P>, Self::Unit>> {
37        let eos_f64 = eos.re();
38        let (temperature, [vapor_density, liquid_density]) =
39            PhaseEquilibrium::pure_p(&eos_f64, self.0, None, Default::default())?;
40
41        let t = temperature.into_reduced();
42        let v1 = 1.0 / liquid_density.to_reduced();
43        let v2 = 1.0 / vapor_density.to_reduced();
44        let p = self.0.into_reduced();
45        let t = Gradient::from(t);
46        let t = t + {
47            let v1 = Gradient::from(v1);
48            let v2 = Gradient::from(v2);
49            let p = Gradient::from(p);
50            let x = E::pure_molefracs();
51
52            let residual_entropy = |v| {
53                let (a, s) = first_derivative(
54                    partial2(
55                        |t, &v, x| eos.lift().residual_helmholtz_energy(t, v, x),
56                        &v,
57                        &x,
58                    ),
59                    t,
60                );
61                (a, -s)
62            };
63            let (a1, s1) = residual_entropy(v1);
64            let (a2, s2) = residual_entropy(v2);
65
66            let ln_rho = (v1 / v2).ln();
67            (p * (v2 - v1) + (a2 - a1 + t * ln_rho)) / (s2 - s1 - ln_rho)
68        };
69        Ok(Temperature::from_reduced(t))
70    }
71}