feos_ad/eos/pcsaft/
pcsaft_binary.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::{jacobian, DualNum, DualVec};
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 binary mixture.
10pub struct PcSaftBinary<const N: usize> {
11    parameters: [[f64; N]; 2],
12    kij: f64,
13}
14
15impl<const N: usize> PcSaftBinary<N> {
16    pub fn new(parameters: [[f64; N]; 2], kij: f64) -> Self {
17        Self { parameters, kij }
18    }
19}
20
21impl<const N: usize> ParametersAD for PcSaftBinary<N> {
22    type Parameters<D: DualNum<f64> + Copy> = ([[D; N]; 2], D);
23
24    fn params<D: DualNum<f64> + Copy>(&self) -> Self::Parameters<D> {
25        (self.parameters.map(|p| p.map(D::from)), D::from(self.kij))
26    }
27
28    fn params_from_inner<D: DualNum<f64> + Copy, D2: DualNum<f64, Inner = D> + Copy>(
29        &(parameters, kij): &Self::Parameters<D>,
30    ) -> Self::Parameters<D2> {
31        (
32            parameters.map(|p| p.map(D2::from_inner)),
33            D2::from_inner(kij),
34        )
35    }
36}
37
38impl<const N: usize> NamedParameters for PcSaftBinary<N> {
39    fn index_parameters_mut<'a, D: DualNum<f64> + Copy>(
40        (_, kij): &'a mut Self::Parameters<D>,
41        index: &str,
42    ) -> &'a mut D {
43        match index {
44            "k_ij" => kij,
45            _ => panic!("{index} is not a valid binary PC-SAFT parameter!"),
46        }
47    }
48}
49
50fn hard_sphere<D: DualNum<f64> + Copy>(
51    [m1, m2]: [D; 2],
52    [x1, x2]: [D; 2],
53    [d1, d2]: [D; 2],
54    density: D,
55) -> (D, [D; 7], D, D, D) {
56    // Packing fractions
57    let zeta = [0, 1, 2, 3].map(|k| (m1 * x1 * d1.powi(k) + m2 * x2 * d2.powi(k)) * FRAC_PI_6);
58    let zeta_23 = zeta[2] / zeta[3];
59    let [zeta0, zeta1, zeta2, zeta3] = zeta.map(|z| z * density);
60    let frac_1mz3 = (-zeta3 + 1.0).recip();
61
62    let eta = zeta3;
63    let eta2 = eta * eta;
64    let eta3 = eta2 * eta;
65    let etas = [
66        D::one(),
67        eta,
68        eta2,
69        eta3,
70        eta2 * eta2,
71        eta2 * eta3,
72        eta3 * eta3,
73    ];
74
75    // hard sphere
76    let hs = (zeta1 * zeta2 * frac_1mz3 * 3.0
77        + zeta2.powi(2) * frac_1mz3.powi(2) * zeta_23
78        + (zeta2 * zeta_23.powi(2) - zeta0) * (-zeta3).ln_1p())
79        / std::f64::consts::FRAC_PI_6;
80
81    (hs, etas, zeta2, zeta3, frac_1mz3)
82}
83
84fn hard_chain<D: DualNum<f64> + Copy>(
85    [m1, m2]: [D; 2],
86    [d1, d2]: [D; 2],
87    [rho1, rho2]: [D; 2],
88    zeta2: D,
89    zeta3: D,
90    frac_1mz3: D,
91) -> D {
92    let c = zeta2 * frac_1mz3 * frac_1mz3;
93    let g_hs = [d1, d2].map(|d| frac_1mz3 + d * c * 1.5 - (d * c).powi(2) * (zeta3 - 1.0) * 0.5);
94    -(rho1 * (m1 - 1.0) * g_hs[0].ln() + rho2 * (m2 - 1.0) * g_hs[1].ln())
95}
96
97#[expect(clippy::too_many_arguments)]
98fn dispersion<D: DualNum<f64> + Copy>(
99    [m1, m2]: [D; 2],
100    [sigma1, sigma2]: [D; 2],
101    [epsilon_k1, epsilon_k2]: [D; 2],
102    kij: D,
103    [x1, x2]: [D; 2],
104    t_inv: D,
105    [rho1, rho2]: [D; 2],
106    etas: [D; 7],
107    frac_1mz3: D,
108) -> (D, [D; 3], [D; 3]) {
109    // binary interactions
110    let m = x1 * m1 + x2 * m2;
111    let epsilon_k11 = epsilon_k1 * t_inv;
112    let epsilon_k12 = (epsilon_k1 * epsilon_k2).sqrt() * t_inv;
113    let epsilon_k22 = epsilon_k2 * t_inv;
114    let sigma11_3 = sigma1.powi(3);
115    let sigma12_3 = ((sigma1 + sigma2) * 0.5).powi(3);
116    let sigma22_3 = sigma2.powi(3);
117    let d11 = rho1 * rho1 * m1 * m1 * epsilon_k11 * sigma11_3;
118    let d12 = rho1 * rho2 * m1 * m2 * epsilon_k12 * sigma12_3 * (-kij + 1.0);
119    let d22 = rho2 * rho2 * m2 * m2 * epsilon_k22 * sigma22_3;
120    let rho1mix = d11 + d12 * 2.0 + d22;
121    let rho2mix = d11 * epsilon_k11 + d12 * epsilon_k12 * (-kij + 1.0) * 2.0 + d22 * epsilon_k22;
122
123    // I1, I2 and C1
124    let mm1 = (m - 1.0) / m;
125    let mm2 = (m - 2.0) / m;
126    let mut i1 = D::zero();
127    let mut i2 = D::zero();
128    for i in 0..7 {
129        i1 += (mm1 * (mm2 * A2[i] + A1[i]) + A0[i]) * etas[i];
130        i2 += (mm1 * (mm2 * B2[i] + B1[i]) + B0[i]) * etas[i];
131    }
132    let eta_m2 = frac_1mz3 * frac_1mz3;
133    let c1 = (m * (etas[1] * 8.0 - etas[2] * 2.0) * eta_m2 * eta_m2 + 1.0
134        - (m - 1.0) * (etas[1] * 20.0 - etas[2] * 27.0 + etas[3] * 12.0 - etas[4] * 2.0)
135            / ((etas[1] - 1.0) * (etas[1] - 2.0)).powi(2))
136    .recip();
137
138    // dispersion
139    let disp = (-rho1mix * i1 * 2.0 - rho2mix * m * c1 * i2) * PI;
140
141    (
142        disp,
143        [sigma11_3, sigma12_3, sigma22_3],
144        [epsilon_k11, epsilon_k12, epsilon_k22],
145    )
146}
147
148#[expect(clippy::too_many_arguments)]
149fn dipoles<D: DualNum<f64> + Copy>(
150    [m1, m2]: [D; 2],
151    [sigma1, sigma2]: [D; 2],
152    [sigma11_3, sigma12_3, sigma22_3]: [D; 3],
153    [epsilon_k11, epsilon_k12, epsilon_k22]: [D; 3],
154    [mu1, mu2]: [D; 2],
155    temperature: D,
156    [rho1, rho2]: [D; 2],
157    etas: [D; 7],
158) -> D {
159    let mu_term1 = mu1 * mu1 / (m1 * temperature * 1.380649e-4) * rho1;
160    let mu_term2 = mu2 * mu2 / (m2 * temperature * 1.380649e-4) * rho2;
161    let sigma111 = sigma11_3;
162    let sigma112 = sigma1 * ((sigma1 + sigma2) * 0.5).powi(2);
163    let sigma122 = sigma2 * ((sigma1 + sigma2) * 0.5).powi(2);
164    let sigma222 = sigma22_3;
165
166    let m11_dipole = if m1.re() > 2.0 { D::from(2.0) } else { m1 };
167    let m22_dipole = if m2.re() > 2.0 { D::from(2.0) } else { m2 };
168    let m12_dipole = (m11_dipole * m22_dipole).sqrt();
169    let [j2_11, j2_12, j2_22] = [
170        (m11_dipole, epsilon_k11),
171        (m12_dipole, epsilon_k12),
172        (m22_dipole, epsilon_k22),
173    ]
174    .map(|(m, e)| {
175        let m1 = (m - 1.0) / m;
176        let m2 = m1 * (m - 2.0) / m;
177        let mut j2 = D::zero();
178        for i in 0..5 {
179            let a = m2 * AD[i][2] + m1 * AD[i][1] + AD[i][0];
180            let b = m2 * BD[i][2] + m1 * BD[i][1] + BD[i][0];
181            j2 += (a + b * e) * etas[i];
182        }
183        j2
184    });
185    let m112_dipole = (m11_dipole * m11_dipole * m22_dipole).cbrt();
186    let m122_dipole = (m11_dipole * m22_dipole * m22_dipole).cbrt();
187    let [j3_111, j3_112, j3_122, j3_222] =
188        [m11_dipole, m112_dipole, m122_dipole, m22_dipole].map(|m| {
189            let m1 = (m - 1.0) / m;
190            let m2 = m1 * (m - 2.0) / m;
191            let mut j3 = D::zero();
192            for i in 0..4 {
193                j3 += (m2 * CD[i][2] + m1 * CD[i][1] + CD[i][0]) * etas[i];
194            }
195            j3
196        });
197
198    let phi2 = (mu_term1 * mu_term1 / sigma11_3 * j2_11
199        + mu_term1 * mu_term2 / sigma12_3 * j2_12 * 2.0
200        + mu_term2 * mu_term2 / sigma22_3 * j2_22)
201        * (-PI);
202    let phi3 = (mu_term1.powi(3) / sigma111 * j3_111
203        + mu_term1.powi(2) * mu_term2 / sigma112 * j3_112 * 3.0
204        + mu_term1 * mu_term2.powi(2) / sigma122 * j3_122 * 3.0
205        + mu_term2.powi(3) / sigma222 * j3_222)
206        * (-PI_SQ_43);
207
208    let mut polar = phi2 * phi2 / (phi2 - phi3);
209    if polar.re().is_nan() {
210        polar = phi2
211    }
212
213    polar
214}
215
216fn association<D: DualNum<f64> + Copy>(
217    assoc_params: [[D; 4]; 2],
218    [sigma11_3, _, sigma22_3]: [D; 3],
219    t_inv: D,
220    [rho1, rho2]: [D; 2],
221    [d1, d2]: [D; 2],
222    zeta2: D,
223    frac_1mz3: D,
224) -> D {
225    let [[kappa_ab1, epsilon_k_ab1, na1, nb1], [kappa_ab2, epsilon_k_ab2, na2, nb2]] = assoc_params;
226
227    let d11 = d1 * 0.5;
228    let d12 = d1 * d2 / (d1 + d2);
229    let d22 = d2 * 0.5;
230    let [k11, k12, k22] = [d11, d12, d22].map(|d| d * zeta2 * frac_1mz3);
231    let s11 = sigma11_3 * kappa_ab1;
232    let mut s12 = (sigma11_3 * sigma22_3 * kappa_ab1 * kappa_ab2).sqrt();
233    if s12.re() == 0.0 {
234        s12 = D::zero();
235    }
236    let s22 = sigma22_3 * kappa_ab2;
237    let e11 = (epsilon_k_ab1 * t_inv).exp() - 1.0;
238    let e12 = ((epsilon_k_ab1 + epsilon_k_ab2) * 0.5 * t_inv).exp() - 1.0;
239    let e22 = (epsilon_k_ab2 * t_inv).exp() - 1.0;
240    let d11 = frac_1mz3 * (k11 * (k11 * 2.0 + 3.0) + 1.0) * s11 * e11;
241    let d12 = frac_1mz3 * (k12 * (k12 * 2.0 + 3.0) + 1.0) * s12 * e12;
242    let d22 = frac_1mz3 * (k22 * (k22 * 2.0 + 3.0) + 1.0) * s22 * e22;
243    let rhoa1 = rho1 * na1;
244    let rhob1 = rho1 * nb1;
245    let rhoa2 = rho2 * na2;
246    let rhob2 = rho2 * nb2;
247
248    let [mut xa1, mut xa2] = [D::from(0.2); 2];
249    for _ in 0..50 {
250        let (g, j) = jacobian(
251            |x| {
252                let [xa1, xa2] = x.data.0[0];
253                let xb1_i =
254                    xa1 * DualVec::from_re(rhoa1 * d11) + xa2 * DualVec::from_re(rhoa2 * d12) + 1.0;
255                let xb2_i =
256                    xa1 * DualVec::from_re(rhoa1 * d12) + xa2 * DualVec::from_re(rhoa2 * d22) + 1.0;
257
258                let f1 = xa1 - 1.0
259                    + xa1 / xb1_i * DualVec::from_re(rhob1 * d11)
260                    + xa1 / xb2_i * DualVec::from_re(rhob2 * d12);
261                let f2 = xa2 - 1.0
262                    + xa2 / xb1_i * DualVec::from_re(rhob1 * d12)
263                    + xa2 / xb2_i * DualVec::from_re(rhob2 * d22);
264
265                SVector::from([f1, f2])
266            },
267            SVector::from([xa1, xa2]),
268        );
269
270        let [g1, g2] = g.data.0[0];
271        let [[j11, j12], [j21, j22]] = j.data.0;
272        let det = j11 * j22 - j12 * j21;
273
274        let delta_xa1 = (j22 * g1 - j12 * g2) / det;
275        let delta_xa2 = (-j21 * g1 + j11 * g2) / det;
276        if delta_xa1.re() < xa1.re() * 0.8 {
277            xa1 -= delta_xa1;
278        } else {
279            xa1 *= 0.2;
280        }
281        if delta_xa2.re() < xa2.re() * 0.8 {
282            xa2 -= delta_xa2;
283        } else {
284            xa2 *= 0.2;
285        }
286
287        if g1.re().abs() < 1e-15 && g2.re().abs() < 1e-15 {
288            break;
289        }
290    }
291
292    let xb1 = (xa1 * rhoa1 * d11 + xa2 * rhoa2 * d12 + 1.0).recip();
293    let xb2 = (xa1 * rhoa1 * d12 + xa2 * rhoa2 * d22 + 1.0).recip();
294    let f = |x: D| x.ln() - x * 0.5 + 0.5;
295
296    rhoa1 * f(xa1) + rhoa2 * f(xa2) + rhob1 * f(xb1) + rhob2 * f(xb2)
297}
298
299#[expect(clippy::too_many_arguments)]
300fn helmholtz_energy_density<D: DualNum<f64> + Copy>(
301    temperature: D,
302    rho: [D; 2],
303    m: [D; 2],
304    sigma: [D; 2],
305    epsilon_k: [D; 2],
306    mu: [D; 2],
307    kij: D,
308    assoc_params: Option<[[D; 4]; 2]>,
309) -> D {
310    // temperature dependent segment diameter
311    let t_inv = temperature.recip();
312    let [sigma1, sigma2] = sigma;
313    let [epsilon_k1, epsilon_k2] = epsilon_k;
314    let d1 = sigma1 * (-(epsilon_k1 * -3. * t_inv).exp() * 0.12 + 1.0);
315    let d2 = sigma2 * (-(epsilon_k2 * -3. * t_inv).exp() * 0.12 + 1.0);
316    let d = [d1, d2];
317
318    // density and composition
319    let [rho1, rho2] = rho;
320    let density = rho1 + rho2;
321    let x = [rho1 / density, rho2 / density];
322
323    // hard sphere
324    let (hs, etas, zeta2, zeta3, frac_1mz3) = hard_sphere(m, x, d, density);
325
326    // hard chain
327    let hc = hard_chain(m, d, rho, zeta2, zeta3, frac_1mz3);
328
329    // dispersion
330    let (disp, sigma_3, epsilon_k_mix) =
331        dispersion(m, sigma, epsilon_k, kij, x, t_inv, rho, etas, frac_1mz3);
332
333    // dipoles
334    let polar = dipoles(m, sigma, sigma_3, epsilon_k_mix, mu, temperature, rho, etas);
335
336    // association
337    let phi = if let Some(p) = assoc_params {
338        let assoc = association(p, sigma_3, t_inv, rho, d, zeta2, frac_1mz3);
339        hs + hc + disp + polar + assoc
340    } else {
341        hs + hc + disp + polar
342    };
343
344    phi * temperature
345}
346
347impl ResidualHelmholtzEnergy<2> for PcSaftBinary<4> {
348    const RESIDUAL: &str = "PC-SAFT (binary)";
349
350    fn compute_max_density(&self, molefracs: &SVector<f64, 2>) -> f64 {
351        let [p1, p2] = self.parameters;
352        let [x1, x2] = molefracs.data.0[0];
353        let [m1, sigma1, ..] = p1;
354        let [m2, sigma2, ..] = p2;
355        MAX_ETA / (FRAC_PI_6 * (m1 * sigma1.powi(3) * x1 + m2 * sigma2.powi(3) * x2))
356    }
357
358    fn residual_helmholtz_energy_density<D: DualNum<f64> + Copy>(
359        &([p1, p2], kij): &Self::Parameters<D>,
360        temperature: D,
361        partial_density: &SVector<D, 2>,
362    ) -> D {
363        let [m1, sigma1, epsilon_k1, mu1] = p1;
364        let [m2, sigma2, epsilon_k2, mu2] = p2;
365        let m = [m1, m2];
366        let sigma = [sigma1, sigma2];
367        let epsilon_k = [epsilon_k1, epsilon_k2];
368        let mu = [mu1, mu2];
369
370        let [rho1, rho2] = partial_density.data.0[0];
371        let rho = [rho1, rho2];
372
373        helmholtz_energy_density(temperature, rho, m, sigma, epsilon_k, mu, kij, None)
374    }
375}
376
377impl ResidualHelmholtzEnergy<2> for PcSaftBinary<8> {
378    const RESIDUAL: &str = "PC-SAFT (binary)";
379
380    fn compute_max_density(&self, molefracs: &SVector<f64, 2>) -> f64 {
381        let [p1, p2] = self.parameters;
382        let [x1, x2] = molefracs.data.0[0];
383        let [m1, sigma1, ..] = p1;
384        let [m2, sigma2, ..] = p2;
385        MAX_ETA / (FRAC_PI_6 * (m1 * sigma1.powi(3) * x1 + m2 * sigma2.powi(3) * x2))
386    }
387
388    fn residual_helmholtz_energy_density<D: DualNum<f64> + Copy>(
389        &([p1, p2], kij): &Self::Parameters<D>,
390        temperature: D,
391        partial_density: &SVector<D, 2>,
392    ) -> D {
393        let [m1, sigma1, epsilon_k1, mu1, kappa_ab1, epsilon_k_ab1, na1, nb1] = p1;
394        let [m2, sigma2, epsilon_k2, mu2, kappa_ab2, epsilon_k_ab2, na2, nb2] = p2;
395        let m = [m1, m2];
396        let sigma = [sigma1, sigma2];
397        let epsilon_k = [epsilon_k1, epsilon_k2];
398        let mu = [mu1, mu2];
399        let assoc_params = Some([
400            [kappa_ab1, epsilon_k_ab1, na1, nb1],
401            [kappa_ab2, epsilon_k_ab2, na2, nb2],
402        ]);
403
404        let [rho1, rho2] = partial_density.data.0[0];
405        let rho = [rho1, rho2];
406
407        helmholtz_energy_density(temperature, rho, m, sigma, epsilon_k, mu, kij, assoc_params)
408    }
409}
410
411#[cfg(test)]
412pub mod test {
413    use super::{PcSaftBinary, ResidualHelmholtzEnergy};
414    use crate::eos::pcsaft::test::pcsaft_binary;
415    use approx::assert_relative_eq;
416    use feos_core::{Contributions::Total, EosResult, ReferenceSystem, State};
417    use nalgebra::SVector;
418    use ndarray::arr1;
419    use quantity::{KELVIN, KILO, METER, MOL};
420
421    #[test]
422    fn test_pcsaft_binary() -> EosResult<()> {
423        let (pcsaft, eos) = pcsaft_binary()?;
424        let pcsaft = (pcsaft.parameters, pcsaft.kij);
425
426        let temperature = 300.0 * KELVIN;
427        let volume = 2.3 * METER * METER * METER;
428        let moles = arr1(&[1.3, 2.5]) * KILO * MOL;
429
430        let state = State::new_nvt(&eos, temperature, volume, &moles)?;
431        let a_feos = state.residual_molar_helmholtz_energy();
432        let mu_feos = state.residual_chemical_potential();
433        let p_feos = state.pressure(Total);
434        let s_feos = state.residual_molar_entropy();
435        let h_feos = state.residual_molar_enthalpy();
436
437        let total_moles = moles.sum();
438        let t = temperature.to_reduced();
439        let v = (volume / total_moles).to_reduced();
440        let x = SVector::from_fn(|i, _| moles.get(i).convert_into(total_moles));
441        let a_ad = PcSaftBinary::residual_molar_helmholtz_energy(&pcsaft, t, v, &x);
442        let mu_ad = PcSaftBinary::residual_chemical_potential(&pcsaft, t, v, &x);
443        let p_ad = PcSaftBinary::pressure(&pcsaft, t, v, &x);
444        let s_ad = PcSaftBinary::residual_molar_entropy(&pcsaft, t, v, &x);
445        let h_ad = PcSaftBinary::residual_molar_enthalpy(&pcsaft, t, v, &x);
446
447        for (s, c) in state.residual_helmholtz_energy_contributions() {
448            println!("{s:20} {}", (c / state.volume).to_reduced());
449        }
450
451        println!("\nMolar Helmholtz energy:\n{}", a_feos.to_reduced(),);
452        println!("{a_ad}");
453        assert_relative_eq!(a_feos.to_reduced(), a_ad, max_relative = 1e-14);
454
455        println!("\nChemical potential:\n{}", mu_feos.get(0).to_reduced());
456        println!("{}", mu_ad[0]);
457        assert_relative_eq!(mu_feos.get(0).to_reduced(), mu_ad[0], max_relative = 1e-14);
458
459        println!("\nPressure:\n{}", p_feos.to_reduced());
460        println!("{p_ad}");
461        assert_relative_eq!(p_feos.to_reduced(), p_ad, max_relative = 1e-14);
462
463        println!("\nMolar entropy:\n{}", s_feos.to_reduced());
464        println!("{s_ad}");
465        assert_relative_eq!(s_feos.to_reduced(), s_ad, max_relative = 1e-14);
466
467        println!("\nMolar enthalpy:\n{}", h_feos.to_reduced());
468        println!("{h_ad}");
469        assert_relative_eq!(h_feos.to_reduced(), h_ad, max_relative = 1e-14);
470
471        Ok(())
472    }
473}
feos_ad/eos/pcsaft/pcsaft_binary.rs

feos_ad/eos/pcsaft/
pcsaft_binary.rs
