feos/pcsaft/eos/
mod.rs

1use super::parameters::PcSaftParameters;
2use crate::association::Association;
3use crate::hard_sphere::{HardSphere, HardSphereProperties};
4use feos_core::parameter::Parameter;
5use feos_core::{
6    Components, EntropyScaling, EosError, EosResult, Molarweight, ReferenceSystem, Residual, State,
7    StateHD,
8};
9use ndarray::Array1;
10use num_dual::DualNum;
11use quantity::*;
12use std::f64::consts::{FRAC_PI_6, PI};
13use std::fmt;
14use std::sync::Arc;
15use typenum::P2;
16
17pub(crate) mod dispersion;
18pub(crate) mod hard_chain;
19pub(crate) mod polar;
20use dispersion::Dispersion;
21use hard_chain::HardChain;
22pub use polar::DQVariants;
23use polar::{Dipole, DipoleQuadrupole, Quadrupole};
24
25/// Customization options for the PC-SAFT equation of state and functional.
26#[derive(Copy, Clone)]
27pub struct PcSaftOptions {
28    pub max_eta: f64,
29    pub max_iter_cross_assoc: usize,
30    pub tol_cross_assoc: f64,
31    pub dq_variant: DQVariants,
32}
33
34impl Default for PcSaftOptions {
35    fn default() -> Self {
36        Self {
37            max_eta: 0.5,
38            max_iter_cross_assoc: 50,
39            tol_cross_assoc: 1e-10,
40            dq_variant: DQVariants::DQ35,
41        }
42    }
43}
44
45/// PC-SAFT equation of state.
46pub struct PcSaft {
47    parameters: Arc<PcSaftParameters>,
48    options: PcSaftOptions,
49    hard_sphere: HardSphere<PcSaftParameters>,
50    hard_chain: Option<HardChain>,
51    dispersion: Dispersion,
52    dipole: Option<Dipole>,
53    quadrupole: Option<Quadrupole>,
54    dipole_quadrupole: Option<DipoleQuadrupole>,
55    association: Option<Association<PcSaftParameters>>,
56}
57
58impl PcSaft {
59    pub fn new(parameters: Arc<PcSaftParameters>) -> Self {
60        Self::with_options(parameters, PcSaftOptions::default())
61    }
62
63    pub fn with_options(parameters: Arc<PcSaftParameters>, options: PcSaftOptions) -> Self {
64        let hard_sphere = HardSphere::new(&parameters);
65        let dispersion = Dispersion {
66            parameters: parameters.clone(),
67        };
68        let hard_chain = if parameters.m.iter().any(|m| (m - 1.0).abs() > 1e-15) {
69            Some(HardChain {
70                parameters: parameters.clone(),
71            })
72        } else {
73            None
74        };
75
76        let dipole = if parameters.ndipole > 0 {
77            Some(Dipole {
78                parameters: parameters.clone(),
79            })
80        } else {
81            None
82        };
83        let quadrupole = if parameters.nquadpole > 0 {
84            Some(Quadrupole {
85                parameters: parameters.clone(),
86            })
87        } else {
88            None
89        };
90        let dipole_quadrupole = if parameters.ndipole > 0 && parameters.nquadpole > 0 {
91            Some(DipoleQuadrupole {
92                parameters: parameters.clone(),
93                variant: options.dq_variant,
94            })
95        } else {
96            None
97        };
98        let association = if !parameters.association.is_empty() {
99            Some(Association::new(
100                &parameters,
101                &parameters.association,
102                options.max_iter_cross_assoc,
103                options.tol_cross_assoc,
104            ))
105        } else {
106            None
107        };
108
109        Self {
110            parameters,
111            options,
112            hard_sphere,
113            hard_chain,
114            dispersion,
115            dipole,
116            quadrupole,
117            dipole_quadrupole,
118            association,
119        }
120    }
121}
122
123impl Components for PcSaft {
124    fn components(&self) -> usize {
125        self.parameters.pure_records.len()
126    }
127
128    fn subset(&self, component_list: &[usize]) -> Self {
129        Self::with_options(
130            Arc::new(self.parameters.subset(component_list)),
131            self.options,
132        )
133    }
134}
135
136impl Residual for PcSaft {
137    fn compute_max_density(&self, moles: &Array1<f64>) -> f64 {
138        self.options.max_eta * moles.sum()
139            / (FRAC_PI_6 * &self.parameters.m * self.parameters.sigma.mapv(|v| v.powi(3)) * moles)
140                .sum()
141    }
142
143    fn residual_helmholtz_energy_contributions<D: DualNum<f64> + Copy>(
144        &self,
145        state: &StateHD<D>,
146    ) -> Vec<(String, D)> {
147        let mut v = Vec::with_capacity(7);
148        let d = self.parameters.hs_diameter(state.temperature);
149
150        v.push((
151            self.hard_sphere.to_string(),
152            self.hard_sphere.helmholtz_energy(state),
153        ));
154        if let Some(hc) = self.hard_chain.as_ref() {
155            v.push((hc.to_string(), hc.helmholtz_energy(state)))
156        }
157        v.push((
158            self.dispersion.to_string(),
159            self.dispersion.helmholtz_energy(state, &d),
160        ));
161        if let Some(dipole) = self.dipole.as_ref() {
162            v.push((dipole.to_string(), dipole.helmholtz_energy(state, &d)))
163        }
164        if let Some(quadrupole) = self.quadrupole.as_ref() {
165            v.push((
166                quadrupole.to_string(),
167                quadrupole.helmholtz_energy(state, &d),
168            ))
169        }
170        if let Some(dipole_quadrupole) = self.dipole_quadrupole.as_ref() {
171            v.push((
172                dipole_quadrupole.to_string(),
173                dipole_quadrupole.helmholtz_energy(state, &d),
174            ))
175        }
176        if let Some(association) = self.association.as_ref() {
177            v.push((
178                association.to_string(),
179                association.helmholtz_energy(state, &d),
180            ))
181        }
182        v
183    }
184}
185
186impl Molarweight for PcSaft {
187    fn molar_weight(&self) -> MolarWeight<Array1<f64>> {
188        self.parameters.molarweight.clone() * GRAM / MOL
189    }
190}
191
192impl fmt::Display for PcSaft {
193    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
194        write!(f, "PC-SAFT")
195    }
196}
197
198fn omega11(t: f64) -> f64 {
199    1.06036 * t.powf(-0.15610)
200        + 0.19300 * (-0.47635 * t).exp()
201        + 1.03587 * (-1.52996 * t).exp()
202        + 1.76474 * (-3.89411 * t).exp()
203}
204
205fn omega22(t: f64) -> f64 {
206    1.16145 * t.powf(-0.14874) + 0.52487 * (-0.77320 * t).exp() + 2.16178 * (-2.43787 * t).exp()
207        - 6.435e-4 * t.powf(0.14874) * (18.0323 * t.powf(-0.76830) - 7.27371).sin()
208}
209
210#[inline]
211fn chapman_enskog_thermal_conductivity(
212    temperature: Temperature,
213    molarweight: MolarWeight,
214    m: f64,
215    sigma: f64,
216    epsilon_k: f64,
217) -> ThermalConductivity {
218    let t = temperature.to_reduced();
219    0.083235 * (t * m / molarweight.convert_to(GRAM / MOL)).sqrt()
220        / sigma.powi(2)
221        / omega22(t / epsilon_k)
222        * WATT
223        / METER
224        / KELVIN
225}
226
227impl EntropyScaling for PcSaft {
228    fn viscosity_reference(
229        &self,
230        temperature: Temperature,
231        _: Volume,
232        moles: &Moles<Array1<f64>>,
233    ) -> EosResult<Viscosity> {
234        let p = &self.parameters;
235        let mw = &p.molarweight;
236        let x = (moles / moles.sum()).into_value();
237        let ce: Array1<_> = (0..self.components())
238            .map(|i| {
239                let tr = (temperature / p.epsilon_k[i] / KELVIN).into_value();
240                5.0 / 16.0 * (mw[i] * GRAM / MOL * KB / NAV * temperature / PI).sqrt()
241                    / omega22(tr)
242                    / (p.sigma[i] * ANGSTROM).powi::<P2>()
243            })
244            .collect();
245        let mut ce_mix = 0.0 * MILLI * PASCAL * SECOND;
246        for i in 0..self.components() {
247            let denom: f64 = (0..self.components())
248                .map(|j| {
249                    x[j] * (1.0
250                        + (ce[i] / ce[j]).into_value().sqrt() * (mw[j] / mw[i]).powf(1.0 / 4.0))
251                    .powi(2)
252                        / (8.0 * (1.0 + mw[i] / mw[j])).sqrt()
253                })
254                .sum();
255            ce_mix += ce[i] * x[i] / denom
256        }
257        Ok(ce_mix)
258    }
259
260    fn viscosity_correlation(&self, s_res: f64, x: &Array1<f64>) -> EosResult<f64> {
261        let coefficients = self
262            .parameters
263            .viscosity
264            .as_ref()
265            .expect("Missing viscosity coefficients.");
266        let m = (x * &self.parameters.m).sum();
267        let s = s_res / m;
268        let pref = (x * &self.parameters.m) / m;
269        let a: f64 = (&coefficients.row(0) * x).sum();
270        let b: f64 = (&coefficients.row(1) * &pref).sum();
271        let c: f64 = (&coefficients.row(2) * &pref).sum();
272        let d: f64 = (&coefficients.row(3) * &pref).sum();
273        Ok(a + b * s + c * s.powi(2) + d * s.powi(3))
274    }
275
276    fn diffusion_reference(
277        &self,
278        temperature: Temperature,
279        volume: Volume,
280        moles: &Moles<Array1<f64>>,
281    ) -> EosResult<Diffusivity> {
282        if self.components() != 1 {
283            return Err(EosError::IncompatibleComponents(self.components(), 1));
284        }
285        let p = &self.parameters;
286        let density = moles.sum() / volume;
287        let res: Array1<_> = (0..self.components())
288            .map(|i| {
289                let tr = (temperature / p.epsilon_k[i] / KELVIN).into_value();
290                3.0 / 8.0 / (p.sigma[i] * ANGSTROM).powi::<P2>() / omega11(tr) / (density * NAV)
291                    * (temperature * RGAS / PI / (p.molarweight[i] * GRAM / MOL) / p.m[i]).sqrt()
292            })
293            .collect();
294        Ok(res[0])
295    }
296
297    fn diffusion_correlation(&self, s_res: f64, x: &Array1<f64>) -> EosResult<f64> {
298        if self.components() != 1 {
299            return Err(EosError::IncompatibleComponents(self.components(), 1));
300        }
301        let coefficients = self
302            .parameters
303            .diffusion
304            .as_ref()
305            .expect("Missing diffusion coefficients.");
306        let m = (x * &self.parameters.m).sum();
307        let s = s_res / m;
308        let pref = (x * &self.parameters.m).mapv(|v| v / m);
309        let a: f64 = (&coefficients.row(0) * x).sum();
310        let b: f64 = (&coefficients.row(1) * &pref).sum();
311        let c: f64 = (&coefficients.row(2) * &pref).sum();
312        let d: f64 = (&coefficients.row(3) * &pref).sum();
313        let e: f64 = (&coefficients.row(4) * &pref).sum();
314        Ok(a + b * s - c * (1.0 - s.exp()) * s.powi(2) - d * s.powi(4) - e * s.powi(8))
315    }
316
317    // Equation 4 of DOI: 10.1021/acs.iecr.9b04289
318    fn thermal_conductivity_reference(
319        &self,
320        temperature: Temperature,
321        volume: Volume,
322        moles: &Moles<Array1<f64>>,
323    ) -> EosResult<ThermalConductivity> {
324        if self.components() != 1 {
325            return Err(EosError::IncompatibleComponents(self.components(), 1));
326        }
327        let p = &self.parameters;
328        let mws = self.molar_weight();
329        let state = State::new_nvt(&Arc::new(Self::new(p.clone())), temperature, volume, moles)?;
330        let res: Array1<_> = (0..self.components())
331            .map(|i| {
332                let tr = (temperature / p.epsilon_k[i] / KELVIN).into_value();
333                let s_res_reduced = state.residual_molar_entropy().to_reduced() / p.m[i];
334                let ref_ce = chapman_enskog_thermal_conductivity(
335                    temperature,
336                    mws.get(i),
337                    p.m[i],
338                    p.sigma[i],
339                    p.epsilon_k[i],
340                );
341                let alpha_visc = (-s_res_reduced / -0.5).exp();
342                let ref_ts = (-0.0167141 * tr / p.m[i] + 0.0470581 * (tr / p.m[i]).powi(2))
343                    * (p.m[i] * p.m[i] * p.sigma[i].powi(3) * p.epsilon_k[i])
344                    * 1e-5
345                    * WATT
346                    / METER
347                    / KELVIN;
348                ref_ce + ref_ts * alpha_visc
349            })
350            .collect();
351        Ok(res[0])
352    }
353
354    fn thermal_conductivity_correlation(&self, s_res: f64, x: &Array1<f64>) -> EosResult<f64> {
355        if self.components() != 1 {
356            return Err(EosError::IncompatibleComponents(self.components(), 1));
357        }
358        let coefficients = self
359            .parameters
360            .thermal_conductivity
361            .as_ref()
362            .expect("Missing thermal conductivity coefficients");
363        let m = (x * &self.parameters.m).sum();
364        let s = s_res / m;
365        let pref = (x * &self.parameters.m).mapv(|v| v / m);
366        let a: f64 = (&coefficients.row(0) * x).sum();
367        let b: f64 = (&coefficients.row(1) * &pref).sum();
368        let c: f64 = (&coefficients.row(2) * &pref).sum();
369        let d: f64 = (&coefficients.row(3) * &pref).sum();
370        Ok(a + b * s + c * (1.0 - s.exp()) + d * s.powi(2))
371    }
372}
373
374#[cfg(test)]
375mod tests {
376    use super::*;
377    use crate::pcsaft::parameters::utils::{
378        butane_parameters, propane_butane_parameters, propane_parameters, water_parameters,
379    };
380    use approx::assert_relative_eq;
381    use feos_core::*;
382    use ndarray::arr1;
383    use quantity::{BAR, KELVIN, METER, MILLI, PASCAL, RGAS, SECOND};
384    use typenum::P3;
385
386    #[test]
387    fn ideal_gas_pressure() {
388        let e = Arc::new(PcSaft::new(propane_parameters()));
389        let t = 200.0 * KELVIN;
390        let v = 1e-3 * METER.powi::<P3>();
391        let n = arr1(&[1.0]) * MOL;
392        let s = State::new_nvt(&e, t, v, &n).unwrap();
393        let p_ig = s.total_moles * RGAS * t / v;
394        assert_relative_eq!(s.pressure(Contributions::IdealGas), p_ig, epsilon = 1e-10);
395        assert_relative_eq!(
396            s.pressure(Contributions::IdealGas) + s.pressure(Contributions::Residual),
397            s.pressure(Contributions::Total),
398            epsilon = 1e-10
399        );
400    }
401
402    #[test]
403    fn ideal_gas_heat_capacity_joback() {
404        let e = Arc::new(PcSaft::new(propane_parameters()));
405        let t = 200.0 * KELVIN;
406        let v = 1e-3 * METER.powi::<P3>();
407        let n = arr1(&[1.0]) * MOL;
408        let s = State::new_nvt(&e, t, v, &n).unwrap();
409        let p_ig = s.total_moles * RGAS * t / v;
410        assert_relative_eq!(s.pressure(Contributions::IdealGas), p_ig, epsilon = 1e-10);
411        assert_relative_eq!(
412            s.pressure(Contributions::IdealGas) + s.pressure(Contributions::Residual),
413            s.pressure(Contributions::Total),
414            epsilon = 1e-10
415        );
416    }
417
418    #[test]
419    fn hard_sphere() {
420        let hs = HardSphere::new(&propane_parameters());
421        let t = 250.0;
422        let v = 1000.0;
423        let n = 1.0;
424        let s = StateHD::new(t, v, arr1(&[n]));
425        let a_rust = hs.helmholtz_energy(&s);
426        assert_relative_eq!(a_rust, 0.410610492598808, epsilon = 1e-10);
427    }
428
429    #[test]
430    fn hard_sphere_mix() {
431        let c1 = HardSphere::new(&propane_parameters());
432        let c2 = HardSphere::new(&butane_parameters());
433        let c12 = HardSphere::new(&propane_butane_parameters());
434        let t = 250.0;
435        let v = 2.5e28;
436        let n = 1.0;
437        let s = StateHD::new(t, v, arr1(&[n]));
438        let a1 = c1.helmholtz_energy(&s);
439        let a2 = c2.helmholtz_energy(&s);
440        let s1m = StateHD::new(t, v, arr1(&[n, 0.0]));
441        let a1m = c12.helmholtz_energy(&s1m);
442        let s2m = StateHD::new(t, v, arr1(&[0.0, n]));
443        let a2m = c12.helmholtz_energy(&s2m);
444        assert_relative_eq!(a1, a1m, epsilon = 1e-14);
445        assert_relative_eq!(a2, a2m, epsilon = 1e-14);
446    }
447
448    #[test]
449    fn association() {
450        let parameters = Arc::new(water_parameters());
451        let assoc = Association::new(&parameters, &parameters.association, 50, 1e-10);
452        let t = 350.0;
453        let v = 41.248289328513216;
454        let n = 1.23;
455        let s = StateHD::new(t, v, arr1(&[n]));
456        let d = parameters.hs_diameter(t);
457        let a_rust = assoc.helmholtz_energy(&s, &d) / n;
458        assert_relative_eq!(a_rust, -4.229878997054543, epsilon = 1e-10);
459    }
460
461    #[test]
462    fn cross_association() {
463        let parameters = Arc::new(water_parameters());
464        let assoc =
465            Association::new_cross_association(&parameters, &parameters.association, 50, 1e-10);
466        let t = 350.0;
467        let v = 41.248289328513216;
468        let n = 1.23;
469        let s = StateHD::new(t, v, arr1(&[n]));
470        let d = parameters.hs_diameter(t);
471        let a_rust = assoc.helmholtz_energy(&s, &d) / n;
472        assert_relative_eq!(a_rust, -4.229878997054543, epsilon = 1e-10);
473    }
474
475    #[test]
476    fn new_tpn() {
477        let e = Arc::new(PcSaft::new(propane_parameters()));
478        let t = 300.0 * KELVIN;
479        let p = BAR;
480        let m = arr1(&[1.0]) * MOL;
481        let s = State::new_npt(&e, t, p, &m, DensityInitialization::None);
482        let p_calc = if let Ok(state) = s {
483            state.pressure(Contributions::Total)
484        } else {
485            0.0 * PASCAL
486        };
487        assert_relative_eq!(p, p_calc, epsilon = 1e-6);
488    }
489
490    #[test]
491    fn vle_pure() {
492        let e = Arc::new(PcSaft::new(propane_parameters()));
493        let t = 300.0 * KELVIN;
494        let vle = PhaseEquilibrium::pure(&e, t, None, Default::default());
495        if let Ok(v) = vle {
496            assert_relative_eq!(
497                v.vapor().pressure(Contributions::Total),
498                v.liquid().pressure(Contributions::Total),
499                epsilon = 1e-6
500            )
501        }
502    }
503
504    #[test]
505    fn critical_point() {
506        let e = Arc::new(PcSaft::new(propane_parameters()));
507        let t = 300.0 * KELVIN;
508        let cp = State::critical_point(&e, None, Some(t), Default::default());
509        if let Ok(v) = cp {
510            assert_relative_eq!(v.temperature, 375.1244078318015 * KELVIN, epsilon = 1e-8)
511        }
512    }
513
514    #[test]
515    fn mix_single() {
516        let e1 = Arc::new(PcSaft::new(propane_parameters()));
517        let e2 = Arc::new(PcSaft::new(butane_parameters()));
518        let e12 = Arc::new(PcSaft::new(propane_butane_parameters()));
519        let t = 300.0 * KELVIN;
520        let v = 0.02456883872966545 * METER.powi::<P3>();
521        let m1 = arr1(&[2.0]) * MOL;
522        let m1m = arr1(&[2.0, 0.0]) * MOL;
523        let m2m = arr1(&[0.0, 2.0]) * MOL;
524        let s1 = State::new_nvt(&e1, t, v, &m1).unwrap();
525        let s2 = State::new_nvt(&e2, t, v, &m1).unwrap();
526        let s1m = State::new_nvt(&e12, t, v, &m1m).unwrap();
527        let s2m = State::new_nvt(&e12, t, v, &m2m).unwrap();
528        assert_relative_eq!(
529            s1.pressure(Contributions::Total),
530            s1m.pressure(Contributions::Total),
531            epsilon = 1e-12
532        );
533        assert_relative_eq!(
534            s2.pressure(Contributions::Total),
535            s2m.pressure(Contributions::Total),
536            epsilon = 1e-12
537        )
538    }
539
540    #[test]
541    fn viscosity() -> EosResult<()> {
542        let e = Arc::new(PcSaft::new(propane_parameters()));
543        let t = 300.0 * KELVIN;
544        let p = BAR;
545        let n = arr1(&[1.0]) * MOL;
546        let s = State::new_npt(&e, t, p, &n, DensityInitialization::None).unwrap();
547        assert_relative_eq!(
548            s.viscosity()?,
549            0.00797 * MILLI * PASCAL * SECOND,
550            epsilon = 1e-5
551        );
552        assert_relative_eq!(
553            s.ln_viscosity_reduced()?,
554            (s.viscosity()? / e.viscosity_reference(s.temperature, s.volume, &s.moles)?)
555                .into_value()
556                .ln(),
557            epsilon = 1e-15
558        );
559        Ok(())
560    }
561
562    #[test]
563    fn diffusion() -> EosResult<()> {
564        let e = Arc::new(PcSaft::new(propane_parameters()));
565        let t = 300.0 * KELVIN;
566        let p = BAR;
567        let n = arr1(&[1.0]) * MOL;
568        let s = State::new_npt(&e, t, p, &n, DensityInitialization::None).unwrap();
569        assert_relative_eq!(
570            s.diffusion()?,
571            0.01505 * (CENTI * METER).powi::<P2>() / SECOND,
572            epsilon = 1e-5
573        );
574        assert_relative_eq!(
575            s.ln_diffusion_reduced()?,
576            (s.diffusion()? / e.diffusion_reference(s.temperature, s.volume, &s.moles)?)
577                .into_value()
578                .ln(),
579            epsilon = 1e-15
580        );
581        Ok(())
582    }
583}
feos/pcsaft/eos/mod.rs

feos/pcsaft/eos/
mod.rs
