feos/epcsaft/eos/
mod.rs

1use super::parameters::{ElectrolytePcSaftParameters, ElectrolytePcSaftPars};
2use crate::association::Association;
3use crate::hard_sphere::{HardSphere, HardSphereProperties};
4use feos_core::{FeosResult, ResidualDyn, Subset};
5use feos_core::{Molarweight, StateHD};
6use nalgebra::DVector;
7use num_dual::DualNum;
8use quantity::*;
9use std::f64::consts::FRAC_PI_6;
10
11pub(crate) mod born;
12pub(crate) mod dispersion;
13pub(crate) mod hard_chain;
14pub(crate) mod ionic;
15pub(crate) mod permittivity;
16use born::Born;
17use dispersion::Dispersion;
18use hard_chain::HardChain;
19use ionic::Ionic;
20
21/// Implemented variants of the ePC-SAFT equation of state.
22#[derive(Copy, Clone, PartialEq)]
23pub enum ElectrolytePcSaftVariants {
24    Advanced,
25    Revised,
26}
27
28/// Customization options for the ePC-SAFT equation of state.
29#[derive(Copy, Clone)]
30pub struct ElectrolytePcSaftOptions {
31    pub max_eta: f64,
32    pub max_iter_cross_assoc: usize,
33    pub tol_cross_assoc: f64,
34    pub epcsaft_variant: ElectrolytePcSaftVariants,
35}
36
37impl Default for ElectrolytePcSaftOptions {
38    fn default() -> Self {
39        Self {
40            max_eta: 0.5,
41            max_iter_cross_assoc: 50,
42            tol_cross_assoc: 1e-10,
43            epcsaft_variant: ElectrolytePcSaftVariants::Advanced,
44        }
45    }
46}
47
48/// electrolyte PC-SAFT (ePC-SAFT) equation of state.
49pub struct ElectrolytePcSaft {
50    pub parameters: ElectrolytePcSaftParameters,
51    pub params: ElectrolytePcSaftPars,
52    pub options: ElectrolytePcSaftOptions,
53    hard_chain: bool,
54    association: Option<Association>,
55    ionic: bool,
56    born: bool,
57}
58
59impl ElectrolytePcSaft {
60    pub fn new(parameters: ElectrolytePcSaftParameters) -> FeosResult<Self> {
61        Self::with_options(parameters, ElectrolytePcSaftOptions::default())
62    }
63
64    pub fn with_options(
65        parameters: ElectrolytePcSaftParameters,
66        options: ElectrolytePcSaftOptions,
67    ) -> FeosResult<Self> {
68        let params = ElectrolytePcSaftPars::new(&parameters)?;
69        let hard_chain = params.m.iter().any(|m| (m - 1.0).abs() > 1e-15);
70
71        let association = (!parameters.association.is_empty())
72            .then(|| Association::new(options.max_iter_cross_assoc, options.tol_cross_assoc));
73
74        let ionic = params.nionic > 0;
75        let born = ionic && matches!(options.epcsaft_variant, ElectrolytePcSaftVariants::Advanced);
76
77        Ok(Self {
78            parameters,
79            params,
80            options,
81            hard_chain,
82            association,
83            ionic,
84            born,
85        })
86    }
87}
88
89impl Subset for ElectrolytePcSaft {
90    fn subset(&self, component_list: &[usize]) -> Self {
91        Self::with_options(self.parameters.subset(component_list), self.options).unwrap()
92    }
93}
94
95impl ResidualDyn for ElectrolytePcSaft {
96    fn components(&self) -> usize {
97        self.parameters.pure.len()
98    }
99
100    fn compute_max_density<D: DualNum<f64> + Copy>(&self, molefracs: &DVector<D>) -> D {
101        let msigma3 = self
102            .params
103            .m
104            .component_mul(&self.params.sigma.map(|v| v.powi(3)));
105        (msigma3.map(D::from).dot(molefracs) * FRAC_PI_6).recip() * self.options.max_eta
106    }
107
108    fn reduced_helmholtz_energy_density_contributions<D: DualNum<f64> + Copy>(
109        &self,
110        state: &StateHD<D>,
111    ) -> Vec<(&'static str, D)> {
112        let mut v = Vec::with_capacity(7);
113        let d = self.params.hs_diameter(state.temperature);
114
115        v.push((
116            "Hard Sphere",
117            HardSphere.helmholtz_energy_density(&self.params, state),
118        ));
119        if self.hard_chain {
120            v.push((
121                "Hard Chain",
122                HardChain.helmholtz_energy_density(&self.params, state),
123            ))
124        }
125        v.push((
126            "Dispersion",
127            Dispersion.helmholtz_energy_density(&self.params, state, &d),
128        ));
129        if let Some(association) = self.association.as_ref() {
130            v.push((
131                "Association",
132                association.helmholtz_energy_density(
133                    &self.params,
134                    &self.parameters.association,
135                    state,
136                    &d,
137                ),
138            ))
139        }
140        if self.ionic {
141            v.push((
142                "Ionic",
143                Ionic.helmholtz_energy_density(
144                    &self.params,
145                    state,
146                    &d,
147                    self.options.epcsaft_variant,
148                ),
149            ))
150        };
151        if self.born {
152            v.push((
153                "Born",
154                Born.helmholtz_energy_density(&self.params, state, &d),
155            ))
156        };
157        v
158    }
159}
160
161impl Molarweight for ElectrolytePcSaft {
162    fn molar_weight(&self) -> MolarWeight<DVector<f64>> {
163        self.parameters.molar_weight.clone()
164    }
165}
166
167#[cfg(test)]
168mod tests {
169    use super::*;
170    use crate::epcsaft::parameters::utils::{
171        butane_parameters, propane_butane_parameters, propane_parameters,
172    };
173    use approx::assert_relative_eq;
174    use feos_core::*;
175    use nalgebra::dvector;
176    use typenum::P3;
177
178    #[test]
179    fn ideal_gas_pressure() {
180        let e = ElectrolytePcSaft::new(propane_parameters()).unwrap();
181        let t = 200.0 * KELVIN;
182        let v = 1e-3 * METER.powi::<P3>();
183        let n = dvector![1.0] * MOL;
184        let s = State::new_nvt(&&e, t, v, &n).unwrap();
185        let p_ig = s.total_moles * RGAS * t / v;
186        assert_relative_eq!(s.pressure(Contributions::IdealGas), p_ig, epsilon = 1e-10);
187        assert_relative_eq!(
188            s.pressure(Contributions::IdealGas) + s.pressure(Contributions::Residual),
189            s.pressure(Contributions::Total),
190            epsilon = 1e-10
191        );
192    }
193
194    #[test]
195    fn ideal_gas_heat_capacity_joback() {
196        let e = ElectrolytePcSaft::new(propane_parameters()).unwrap();
197        let t = 200.0 * KELVIN;
198        let v = 1e-3 * METER.powi::<P3>();
199        let n = dvector![1.0] * MOL;
200        let s = State::new_nvt(&&e, t, v, &n).unwrap();
201        let p_ig = s.total_moles * RGAS * t / v;
202        assert_relative_eq!(s.pressure(Contributions::IdealGas), p_ig, epsilon = 1e-10);
203        assert_relative_eq!(
204            s.pressure(Contributions::IdealGas) + s.pressure(Contributions::Residual),
205            s.pressure(Contributions::Total),
206            epsilon = 1e-10
207        );
208    }
209
210    #[test]
211    fn hard_sphere() {
212        let p = ElectrolytePcSaftPars::new(&propane_parameters()).unwrap();
213        let t = 250.0;
214        let v = 1000.0;
215        let n = 1.0;
216        let s = StateHD::new(t, v, &dvector![n]);
217        let a_rust = HardSphere.helmholtz_energy_density(&p, &s) * v;
218        assert_relative_eq!(a_rust, 0.410610492598808, epsilon = 1e-10);
219    }
220
221    #[test]
222    fn hard_sphere_mix() {
223        let p1 = ElectrolytePcSaftPars::new(&propane_parameters()).unwrap();
224        let p2 = ElectrolytePcSaftPars::new(&butane_parameters()).unwrap();
225        let p12 = ElectrolytePcSaftPars::new(&propane_butane_parameters()).unwrap();
226        let t = 250.0;
227        let v = 2.5e28;
228        let n = 1.0;
229        let s = StateHD::new(t, v, &dvector![n]);
230        let a1 = HardSphere.helmholtz_energy_density(&p1, &s);
231        let a2 = HardSphere.helmholtz_energy_density(&p2, &s);
232        let s1m = StateHD::new(t, v, &dvector![n, 0.0]);
233        let a1m = HardSphere.helmholtz_energy_density(&p12, &s1m);
234        let s2m = StateHD::new(t, v, &dvector![0.0, n]);
235        let a2m = HardSphere.helmholtz_energy_density(&p12, &s2m);
236        assert_relative_eq!(a1, a1m, epsilon = 1e-14);
237        assert_relative_eq!(a2, a2m, epsilon = 1e-14);
238    }
239
240    #[test]
241    fn new_tpn() {
242        let e = ElectrolytePcSaft::new(propane_parameters()).unwrap();
243        let t = 300.0 * KELVIN;
244        let p = BAR;
245        let m = dvector![1.0] * MOL;
246        let s = State::new_npt(&&e, t, p, &m, None);
247        let p_calc = if let Ok(state) = s {
248            state.pressure(Contributions::Total)
249        } else {
250            0.0 * PASCAL
251        };
252        assert_relative_eq!(p, p_calc, epsilon = 1e-6);
253    }
254
255    #[test]
256    fn vle_pure() {
257        let e = ElectrolytePcSaft::new(propane_parameters()).unwrap();
258        let t = 300.0 * KELVIN;
259        let vle = PhaseEquilibrium::pure(&&e, t, None, Default::default());
260        if let Ok(v) = vle {
261            assert_relative_eq!(
262                v.vapor().pressure(Contributions::Total),
263                v.liquid().pressure(Contributions::Total),
264                epsilon = 1e-6
265            )
266        }
267    }
268
269    #[test]
270    fn critical_point() {
271        let e = ElectrolytePcSaft::new(propane_parameters()).unwrap();
272        let t = 300.0 * KELVIN;
273        let cp = State::critical_point(&&e, None, Some(t), None, Default::default());
274        if let Ok(v) = cp {
275            assert_relative_eq!(v.temperature, 375.1244078318015 * KELVIN, epsilon = 1e-8)
276        }
277    }
278
279    #[test]
280    fn mix_single() {
281        let e1 = ElectrolytePcSaft::new(propane_parameters()).unwrap();
282        let e2 = ElectrolytePcSaft::new(butane_parameters()).unwrap();
283        let e12 = ElectrolytePcSaft::new(propane_butane_parameters()).unwrap();
284        let t = 300.0 * KELVIN;
285        let v = 0.02456883872966545 * METER.powi::<P3>();
286        let m1 = dvector![2.0] * MOL;
287        let m1m = dvector![2.0, 0.0] * MOL;
288        let m2m = dvector![0.0, 2.0] * MOL;
289        let s1 = State::new_nvt(&&e1, t, v, &m1).unwrap();
290        let s2 = State::new_nvt(&&e2, t, v, &m1).unwrap();
291        let s1m = State::new_nvt(&&e12, t, v, &m1m).unwrap();
292        let s2m = State::new_nvt(&&e12, t, v, &m2m).unwrap();
293        assert_relative_eq!(
294            s1.pressure(Contributions::Total),
295            s1m.pressure(Contributions::Total),
296            epsilon = 1e-12
297        );
298        assert_relative_eq!(
299            s2.pressure(Contributions::Total),
300            s2m.pressure(Contributions::Total),
301            epsilon = 1e-12
302        );
303        assert_relative_eq!(
304            s2.pressure(Contributions::Total),
305            s2m.pressure(Contributions::Total),
306            epsilon = 1e-12
307        )
308    }
309}
feos/epcsaft/eos/mod.rs

feos/epcsaft/eos/
mod.rs
