feos_core/equation_of_state/
residual.rs

1use super::Components;
2use crate::{EosError, EosResult, ReferenceSystem, StateHD};
3use ndarray::prelude::*;
4use ndarray::ScalarOperand;
5use num_dual::*;
6use num_traits::{One, Zero};
7use quantity::*;
8use std::ops::Div;
9use typenum::Quot;
10
11/// Molar weight of all components.
12///
13/// Enables calculation of (mass) specific properties.
14pub trait Molarweight {
15    fn molar_weight(&self) -> MolarWeight<Array1<f64>>;
16}
17
18/// A residual Helmholtz energy model.
19pub trait Residual: Components + Send + Sync {
20    /// Return the maximum density in Angstrom^-3.
21    ///
22    /// This value is used as an estimate for a liquid phase for phase
23    /// equilibria and other iterations. It is not explicitly meant to
24    /// be a mathematical limit for the density (if those exist in the
25    /// equation of state anyways).
26    fn compute_max_density(&self, moles: &Array1<f64>) -> f64;
27
28    /// Evaluate the reduced Helmholtz energy of each individual contribution
29    /// and return them together with a string representation of the contribution.
30    fn residual_helmholtz_energy_contributions<D: DualNum<f64> + Copy + ScalarOperand>(
31        &self,
32        state: &StateHD<D>,
33    ) -> Vec<(String, D)>;
34
35    /// Evaluate the residual reduced Helmholtz energy $\beta A^\mathrm{res}$.
36    fn residual_helmholtz_energy<D: DualNum<f64> + Copy + ScalarOperand>(
37        &self,
38        state: &StateHD<D>,
39    ) -> D {
40        self.residual_helmholtz_energy_contributions(state)
41            .iter()
42            .fold(D::zero(), |acc, (_, a)| acc + a)
43    }
44
45    /// Check if the provided optional mole number is consistent with the
46    /// equation of state.
47    ///
48    /// In general, the number of elements in `moles` needs to match the number
49    /// of components of the equation of state. For a pure component, however,
50    /// no moles need to be provided. In that case, it is set to the constant
51    /// reference value.
52    fn validate_moles(&self, moles: Option<&Moles<Array1<f64>>>) -> EosResult<Moles<Array1<f64>>> {
53        let l = moles.map_or(1, |m| m.len());
54        if self.components() == l {
55            match moles {
56                Some(m) => Ok(m.to_owned()),
57                None => Ok(Moles::from_reduced(Array::ones(1))),
58            }
59        } else {
60            Err(EosError::IncompatibleComponents(self.components(), l))
61        }
62    }
63
64    /// Calculate the maximum density.
65    ///
66    /// This value is used as an estimate for a liquid phase for phase
67    /// equilibria and other iterations. It is not explicitly meant to
68    /// be a mathematical limit for the density (if those exist in the
69    /// equation of state anyways).
70    fn max_density(&self, moles: Option<&Moles<Array1<f64>>>) -> EosResult<Density> {
71        let mr = self.validate_moles(moles)?.to_reduced();
72        Ok(Density::from_reduced(self.compute_max_density(&mr)))
73    }
74
75    /// Calculate the second virial coefficient $B(T)$
76    fn second_virial_coefficient(
77        &self,
78        temperature: Temperature,
79        moles: Option<&Moles<Array1<f64>>>,
80    ) -> EosResult<Quot<f64, Density>> {
81        let mr = self.validate_moles(moles)?;
82        let x = (&mr / mr.sum()).into_value();
83        let mut rho = HyperDual64::zero();
84        rho.eps1 = 1.0;
85        rho.eps2 = 1.0;
86        let t = HyperDual64::from(temperature.to_reduced());
87        let s = StateHD::new_virial(t, rho, x);
88        Ok(Quantity::from_reduced(
89            self.residual_helmholtz_energy(&s).eps1eps2 * 0.5,
90        ))
91    }
92
93    /// Calculate the third virial coefficient $C(T)$
94    #[allow(clippy::type_complexity)]
95    fn third_virial_coefficient(
96        &self,
97        temperature: Temperature,
98        moles: Option<&Moles<Array1<f64>>>,
99    ) -> EosResult<<<f64 as Div<Density>>::Output as Div<Density>>::Output> {
100        let mr = self.validate_moles(moles)?;
101        let x = (&mr / mr.sum()).into_value();
102        let rho = Dual3_64::zero().derivative();
103        let t = Dual3_64::from(temperature.to_reduced());
104        let s = StateHD::new_virial(t, rho, x);
105        Ok(Quantity::from_reduced(
106            self.residual_helmholtz_energy(&s).v3 / 3.0,
107        ))
108    }
109
110    /// Calculate the temperature derivative of the second virial coefficient $B'(T)$
111    #[allow(clippy::type_complexity)]
112    fn second_virial_coefficient_temperature_derivative(
113        &self,
114        temperature: Temperature,
115        moles: Option<&Moles<Array1<f64>>>,
116    ) -> EosResult<<<f64 as Div<Density>>::Output as Div<Temperature>>::Output> {
117        let mr = self.validate_moles(moles)?;
118        let x = (&mr / mr.sum()).into_value();
119        let mut rho = HyperDual::zero();
120        rho.eps1 = Dual64::one();
121        rho.eps2 = Dual64::one();
122        let t = HyperDual::from_re(Dual64::from(temperature.to_reduced()).derivative());
123        let s = StateHD::new_virial(t, rho, x);
124        Ok(Quantity::from_reduced(
125            self.residual_helmholtz_energy(&s).eps1eps2.eps * 0.5,
126        ))
127    }
128
129    /// Calculate the temperature derivative of the third virial coefficient $C'(T)$
130    #[allow(clippy::type_complexity)]
131    fn third_virial_coefficient_temperature_derivative(
132        &self,
133        temperature: Temperature,
134        moles: Option<&Moles<Array1<f64>>>,
135    ) -> EosResult<
136        <<<f64 as Div<Density>>::Output as Div<Density>>::Output as Div<Temperature>>::Output,
137    > {
138        let mr = self.validate_moles(moles)?;
139        let x = (&mr / mr.sum()).into_value();
140        let rho = Dual3::zero().derivative();
141        let t = Dual3::from_re(Dual64::from(temperature.to_reduced()).derivative());
142        let s = StateHD::new_virial(t, rho, x);
143        Ok(Quantity::from_reduced(
144            self.residual_helmholtz_energy(&s).v3.eps / 3.0,
145        ))
146    }
147}
148
149/// Reference values and residual entropy correlations for entropy scaling.
150pub trait EntropyScaling {
151    fn viscosity_reference(
152        &self,
153        temperature: Temperature,
154        volume: Volume,
155        moles: &Moles<Array1<f64>>,
156    ) -> EosResult<Viscosity>;
157    fn viscosity_correlation(&self, s_res: f64, x: &Array1<f64>) -> EosResult<f64>;
158    fn diffusion_reference(
159        &self,
160        temperature: Temperature,
161        volume: Volume,
162        moles: &Moles<Array1<f64>>,
163    ) -> EosResult<Diffusivity>;
164    fn diffusion_correlation(&self, s_res: f64, x: &Array1<f64>) -> EosResult<f64>;
165    fn thermal_conductivity_reference(
166        &self,
167        temperature: Temperature,
168        volume: Volume,
169        moles: &Moles<Array1<f64>>,
170    ) -> EosResult<ThermalConductivity>;
171    fn thermal_conductivity_correlation(&self, s_res: f64, x: &Array1<f64>) -> EosResult<f64>;
172}
173
174/// Dummy implementation for [EquationOfState](super::EquationOfState)s that only contain an ideal gas contribution.
175pub struct NoResidual(pub usize);
176
177impl Components for NoResidual {
178    fn components(&self) -> usize {
179        self.0
180    }
181
182    fn subset(&self, component_list: &[usize]) -> Self {
183        Self(component_list.len())
184    }
185}
186
187impl Residual for NoResidual {
188    fn compute_max_density(&self, _: &Array1<f64>) -> f64 {
189        1.0
190    }
191
192    fn residual_helmholtz_energy_contributions<D: DualNum<f64> + Copy>(
193        &self,
194        _: &StateHD<D>,
195    ) -> Vec<(String, D)> {
196        vec![]
197    }
198}
feos_core/equation_of_state/residual.rs

feos_core/equation_of_state/
residual.rs
