feos_pets/eos/
hard_sphere.rs

1use crate::parameters::PetsParameters;
2use feos_core::{HelmholtzEnergyDual, StateHD};
3use ndarray::*;
4use num_dual::DualNum;
5use std::fmt;
6use std::rc::Rc;
7
8impl PetsParameters {
9    pub fn hs_diameter<D: DualNum<f64>>(&self, temperature: D) -> Array1<D> {
10        let ti = temperature.recip() * -3.052785558;
11        Array::from_shape_fn(self.sigma.len(), |i| {
12            -((ti * self.epsilon_k[i]).exp() * 0.127112544 - 1.0) * self.sigma[i]
13        })
14    }
15}
16
17#[derive(Debug, Clone)]
18pub struct HardSphere {
19    pub parameters: Rc<PetsParameters>,
20}
21
22impl<D: DualNum<f64>> HelmholtzEnergyDual<D> for HardSphere {
23    fn helmholtz_energy(&self, state: &StateHD<D>) -> D {
24        let d = self.parameters.hs_diameter(state.temperature);
25        let zeta = zeta(&state.partial_density, &d);
26        let frac_1mz3 = -(zeta[3] - 1.0).recip();
27        let zeta_23 = zeta_23(&state.molefracs, &d);
28
29        state.volume * 6.0 / std::f64::consts::PI
30            * (zeta[1] * zeta[2] * frac_1mz3 * 3.0
31                + zeta[2].powi(2) * frac_1mz3.powi(2) * zeta_23
32                + (zeta[2] * zeta_23.powi(2) - zeta[0]) * (zeta[3] * (-1.0)).ln_1p())
33    }
34}
35
36impl fmt::Display for HardSphere {
37    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
38        write!(f, "Hard Sphere")
39    }
40}
41
42pub fn zeta<D: DualNum<f64>>(partial_density: &Array1<D>, diameter: &Array1<D>) -> [D; 4] {
43    let mut zeta: [D; 4] = [D::zero(), D::zero(), D::zero(), D::zero()];
44    for i in 0..diameter.len() {
45        for (k, z) in zeta.iter_mut().enumerate() {
46            *z += partial_density[i] * diameter[i].powi(k as i32) * (std::f64::consts::PI / 6.0);
47        }
48    }
49    zeta
50}
51
52pub fn zeta_23<D: DualNum<f64>>(molefracs: &Array1<D>, diameter: &Array1<D>) -> D {
53    let mut zeta: [D; 2] = [D::zero(), D::zero()];
54    for i in 0..diameter.len() {
55        for (k, z) in zeta.iter_mut().enumerate() {
56            *z += molefracs[i] * diameter[i].powi((k + 2) as i32);
57        }
58    }
59    zeta[0] / zeta[1]
60}
61
62#[cfg(test)]
63mod tests {
64    use super::*;
65    use crate::parameters::utils::{
66        argon_krypton_parameters, argon_parameters, krypton_parameters,
67    };
68    use approx::assert_relative_eq;
69    use ndarray::arr1;
70
71    #[test]
72    fn mix() {
73        let c1 = HardSphere {
74            parameters: argon_parameters(),
75        };
76        let c2 = HardSphere {
77            parameters: krypton_parameters(),
78        };
79        let c12 = HardSphere {
80            parameters: argon_krypton_parameters(),
81        };
82        let t = 250.0;
83        let v = 2.5e28;
84        let n = 1.0;
85        let s = StateHD::new(t, v, arr1(&[n]));
86        let a1 = c1.helmholtz_energy(&s);
87        let a2 = c2.helmholtz_energy(&s);
88        let s1m = StateHD::new(t, v, arr1(&[n, 0.0]));
89        let a1m = c12.helmholtz_energy(&s1m);
90        let s2m = StateHD::new(t, v, arr1(&[0.0, n]));
91        let a2m = c12.helmholtz_energy(&s2m);
92        assert_relative_eq!(a1, a1m, epsilon = 1e-14);
93        assert_relative_eq!(a2, a2m, epsilon = 1e-14);
94    }
95}
feos_pets/eos/hard_sphere.rs

feos_pets/eos/
hard_sphere.rs
