feos_ad/eos/pcsaft/
pcsaft_pure.rs

1use super::{A0, A1, A2, AD, B0, B1, B2, BD, CD, MAX_ETA};
2use crate::{NamedParameters, ParametersAD, ResidualHelmholtzEnergy};
3use nalgebra::SVector;
4use num_dual::DualNum;
5use std::f64::consts::{FRAC_PI_6, PI};
6
7const PI_SQ_43: f64 = 4.0 / 3.0 * PI * PI;
8
9/// Optimized implementation of PC-SAFT for a single component.
10#[derive(Clone, Copy)]
11pub struct PcSaftPure<const N: usize>(pub [f64; N]);
12
13fn helmholtz_energy_density_non_assoc<D: DualNum<f64> + Copy>(
14    m: D,
15    sigma: D,
16    epsilon_k: D,
17    mu: D,
18    temperature: D,
19    density: D,
20) -> (D, [D; 2]) {
21    // temperature dependent segment diameter
22    let diameter = sigma * (-(epsilon_k * (-3.) / temperature).exp() * 0.12 + 1.0);
23
24    let eta = m * density * diameter.powi(3) * FRAC_PI_6;
25    let eta2 = eta * eta;
26    let eta3 = eta2 * eta;
27    let eta_m1 = (-eta + 1.0).recip();
28    let eta_m2 = eta_m1 * eta_m1;
29    let etas = [
30        D::one(),
31        eta,
32        eta2,
33        eta3,
34        eta2 * eta2,
35        eta2 * eta3,
36        eta3 * eta3,
37    ];
38
39    // hard sphere
40    let hs = m * density * (eta * 4.0 - eta2 * 3.0) * eta_m2;
41
42    // hard chain
43    let g = (-eta * 0.5 + 1.0) * eta_m1 * eta_m2;
44    let hc = -density * (m - 1.0) * g.ln();
45
46    // dispersion
47    let e = epsilon_k / temperature;
48    let s3 = sigma.powi(3);
49    let mut i1 = D::zero();
50    let mut i2 = D::zero();
51    let m1 = (m - 1.0) / m;
52    let m2 = (m - 2.0) / m;
53    for i in 0..7 {
54        i1 += (m1 * (m2 * A2[i] + A1[i]) + A0[i]) * etas[i];
55        i2 += (m1 * (m2 * B2[i] + B1[i]) + B0[i]) * etas[i];
56    }
57    let c1 = (m * (eta * 8.0 - eta2 * 2.0) * eta_m2 * eta_m2 + 1.0
58        - (m - 1.0) * (eta * 20.0 - eta2 * 27.0 + eta2 * eta * 12.0 - eta2 * eta2 * 2.0)
59            / ((eta - 1.0) * (eta - 2.0)).powi(2))
60    .recip();
61    let i = i1 * 2.0 + c1 * i2 * m * e;
62    let disp = -density * density * m.powi(2) * e * s3 * i * PI;
63
64    // dipoles
65    let mu2 = mu.powi(2) / (m * temperature * 1.380649e-4);
66    let m_dipole = if m.re() > 2.0 { D::from(2.0) } else { m };
67    let m1 = (m_dipole - 1.0) / m_dipole;
68    let m2 = m1 * (m_dipole - 2.0) / m_dipole;
69    let mut j1 = D::zero();
70    let mut j2 = D::zero();
71    for i in 0..5 {
72        let a = m2 * AD[i][2] + m1 * AD[i][1] + AD[i][0];
73        let b = m2 * BD[i][2] + m1 * BD[i][1] + BD[i][0];
74        j1 += (a + b * e) * etas[i];
75        if i < 4 {
76            j2 += (m2 * CD[i][2] + m1 * CD[i][1] + CD[i][0]) * etas[i];
77        }
78    }
79
80    // mu is factored out of these expressions to deal with the case where mu=0
81    let phi2 = -density * density * j1 / s3 * PI;
82    let phi3 = -density * density * density * j2 / s3 * PI_SQ_43;
83    let dipole = phi2 * phi2 * mu2 * mu2 / (phi2 - phi3 * mu2);
84
85    ((hs + hc + disp + dipole) * temperature, [eta, eta_m1])
86}
87
88fn helmholtz_energy_density<D: DualNum<f64> + Copy>(
89    parameters: &[D; 8],
90    temperature: D,
91    density: D,
92) -> D {
93    let [m, sigma, epsilon_k, mu, kappa_ab, epsilon_k_ab, na, nb] = *parameters;
94    let (non_assoc, [eta, eta_m1]) =
95        helmholtz_energy_density_non_assoc(m, sigma, epsilon_k, mu, temperature, density);
96
97    // association
98    let delta_assoc = ((epsilon_k_ab / temperature).exp() - 1.0) * sigma.powi(3) * kappa_ab;
99    let k = eta * eta_m1;
100    let delta = (k * (k * 0.5 + 1.5) + 1.0) * eta_m1 * delta_assoc;
101    let rhoa = na * density;
102    let rhob = nb * density;
103    let aux = (rhoa - rhob) * delta + 1.0;
104    let sqrt = (aux * aux + rhob * delta * 4.0).sqrt();
105    let xa = (sqrt + 1.0 + (rhob - rhoa) * delta).recip() * 2.0;
106    let xb = (sqrt + 1.0 - (rhob - rhoa) * delta).recip() * 2.0;
107    let assoc = rhoa * (xa.ln() - xa * 0.5 + 0.5) + rhob * (xb.ln() - xb * 0.5 + 0.5);
108
109    non_assoc + assoc * temperature
110}
111
112impl<const N: usize> ParametersAD for PcSaftPure<N> {
113    type Parameters<D: DualNum<f64> + Copy> = [D; N];
114
115    fn params<D: DualNum<f64> + Copy>(&self) -> Self::Parameters<D> {
116        self.0.map(D::from)
117    }
118
119    fn params_from_inner<D: DualNum<f64> + Copy, D2: DualNum<f64, Inner = D> + Copy>(
120        parameters: &Self::Parameters<D>,
121    ) -> Self::Parameters<D2> {
122        parameters.map(D2::from_inner)
123    }
124}
125
126impl ResidualHelmholtzEnergy<1> for PcSaftPure<8> {
127    const RESIDUAL: &str = "PC-SAFT (pure)";
128
129    fn compute_max_density(&self, _: &SVector<f64, 1>) -> f64 {
130        let m = self.0[0];
131        let sigma = self.0[1];
132        MAX_ETA / (FRAC_PI_6 * m * sigma.powi(3))
133    }
134
135    fn residual_helmholtz_energy_density<D: DualNum<f64> + Copy>(
136        parameters: &Self::Parameters<D>,
137        temperature: D,
138        partial_density: &SVector<D, 1>,
139    ) -> D {
140        let density = partial_density.data.0[0][0];
141        helmholtz_energy_density(parameters, temperature, density)
142    }
143}
144
145impl ResidualHelmholtzEnergy<1> for PcSaftPure<4> {
146    const RESIDUAL: &str = "PC-SAFT (pure)";
147
148    fn compute_max_density(&self, _: &SVector<f64, 1>) -> f64 {
149        let m = self.0[0];
150        let sigma = self.0[1];
151        MAX_ETA / (FRAC_PI_6 * m * sigma.powi(3))
152    }
153
154    fn residual_helmholtz_energy_density<D: DualNum<f64> + Copy>(
155        parameters: &Self::Parameters<D>,
156        temperature: D,
157        partial_density: &SVector<D, 1>,
158    ) -> D {
159        let density = partial_density.data.0[0][0];
160        let [m, sigma, epsilon_k, mu] = *parameters;
161        helmholtz_energy_density_non_assoc(m, sigma, epsilon_k, mu, temperature, density).0
162    }
163}
164
165impl<const N: usize> NamedParameters for PcSaftPure<N> {
166    fn index_parameters_mut<'a, D: DualNum<f64> + Copy>(
167        parameters: &'a mut [D; N],
168        index: &str,
169    ) -> &'a mut D {
170        match index {
171            "m" => &mut parameters[0],
172            "sigma" => &mut parameters[1],
173            "epsilon_k" => &mut parameters[2],
174            "mu" => &mut parameters[3],
175            "kappa_ab" => &mut parameters[4],
176            "epsilon_k_ab" => &mut parameters[5],
177            "na" => &mut parameters[6],
178            "nb" => &mut parameters[7],
179            _ => panic!("{index} is not a valid PC-SAFT parameter!"),
180        }
181    }
182}
183
184#[cfg(test)]
185pub mod test {
186    use super::{PcSaftPure, ResidualHelmholtzEnergy};
187    use crate::eos::pcsaft::test::pcsaft;
188    use approx::assert_relative_eq;
189    use feos_core::{Contributions::Total, EosResult, ReferenceSystem, State};
190    use nalgebra::SVector;
191    use ndarray::arr1;
192    use quantity::{KELVIN, KILO, METER, MOL};
193
194    #[test]
195    fn test_pcsaft_pure() -> EosResult<()> {
196        let (pcsaft, eos) = pcsaft()?;
197        let pcsaft = pcsaft.0;
198
199        let temperature = 300.0 * KELVIN;
200        let volume = 2.3 * METER * METER * METER;
201        let moles = arr1(&[1.3]) * KILO * MOL;
202
203        let state = State::new_nvt(&eos, temperature, volume, &moles)?;
204        let a_feos = state.residual_molar_helmholtz_energy();
205        let mu_feos = state.residual_chemical_potential();
206        let p_feos = state.pressure(Total);
207        let s_feos = state.residual_molar_entropy();
208        let h_feos = state.residual_molar_enthalpy();
209
210        let total_moles = moles.sum();
211        let t = temperature.to_reduced();
212        let v = (volume / total_moles).to_reduced();
213        let x = SVector::from_fn(|i, _| moles.get(i).convert_into(total_moles));
214        let a_ad = PcSaftPure::residual_molar_helmholtz_energy(&pcsaft, t, v, &x);
215        let mu_ad = PcSaftPure::residual_chemical_potential(&pcsaft, t, v, &x);
216        let p_ad = PcSaftPure::pressure(&pcsaft, t, v, &x);
217        let s_ad = PcSaftPure::residual_molar_entropy(&pcsaft, t, v, &x);
218        let h_ad = PcSaftPure::residual_molar_enthalpy(&pcsaft, t, v, &x);
219
220        println!("\nMolar Helmholtz energy:\n{}", a_feos.to_reduced(),);
221        println!("{a_ad}");
222        assert_relative_eq!(a_feos.to_reduced(), a_ad, max_relative = 1e-14);
223
224        println!("\nChemical potential:\n{}", mu_feos.get(0).to_reduced());
225        println!("{}", mu_ad[0]);
226        assert_relative_eq!(mu_feos.get(0).to_reduced(), mu_ad[0], max_relative = 1e-14);
227
228        println!("\nPressure:\n{}", p_feos.to_reduced());
229        println!("{p_ad}");
230        assert_relative_eq!(p_feos.to_reduced(), p_ad, max_relative = 1e-14);
231
232        println!("\nMolar entropy:\n{}", s_feos.to_reduced());
233        println!("{s_ad}");
234        assert_relative_eq!(s_feos.to_reduced(), s_ad, max_relative = 1e-14);
235
236        println!("\nMolar enthalpy:\n{}", h_feos.to_reduced());
237        println!("{h_ad}");
238        assert_relative_eq!(h_feos.to_reduced(), h_ad, max_relative = 1e-14);
239
240        Ok(())
241    }
242}
feos_ad/eos/pcsaft/pcsaft_pure.rs

feos_ad/eos/pcsaft/
pcsaft_pure.rs
