feos_pcsaft/eos/
polar.rs

1use crate::parameters::PcSaftParameters;
2use feos_core::{HelmholtzEnergyDual, StateHD};
3use ndarray::prelude::*;
4use num_dual::DualNum;
5use std::f64::consts::{FRAC_PI_3, PI};
6use std::fmt;
7use std::rc::Rc;
8
9pub const ALPHA: f64 = 1.1937350;
10
11// Dipole parameters
12pub const AD: [[f64; 3]; 5] = [
13    [0.30435038064, 0.95346405973, -1.16100802773],
14    [-0.13585877707, -1.83963831920, 4.52586067320],
15    [1.44933285154, 2.01311801180, 0.97512223853],
16    [0.35569769252, -7.37249576667, -12.2810377713],
17    [-2.06533084541, 8.23741345333, 5.93975747420],
18];
19
20pub const BD: [[f64; 3]; 5] = [
21    [0.21879385627, -0.58731641193, 3.48695755800],
22    [-1.18964307357, 1.24891317047, -14.9159739347],
23    [1.16268885692, -0.50852797392, 15.3720218600],
24    [0.0; 3],
25    [0.0; 3],
26];
27
28pub const CD: [[f64; 3]; 4] = [
29    [-0.06467735252, -0.95208758351, -0.62609792333],
30    [0.19758818347, 2.99242575222, 1.29246858189],
31    [-0.80875619458, -2.38026356489, 1.65427830900],
32    [0.69028490492, -0.27012609786, -3.43967436378],
33];
34
35// Quadrupole parameters
36pub const AQ: [[f64; 3]; 5] = [
37    [1.237830788, 1.285410878, 1.794295401],
38    [2.435503144, -11.46561451, 0.769510293],
39    [1.633090469, 22.08689285, 7.264792255],
40    [-1.611815241, 7.46913832, 94.48669892],
41    [6.977118504, -17.19777208, -77.1484579],
42];
43
44pub const BQ: [[f64; 3]; 5] = [
45    [0.454271755, -0.813734006, 6.868267516],
46    [-4.501626435, 10.06402986, -5.173223765],
47    [3.585886783, -10.87663092, -17.2402066],
48    [0.0; 3],
49    [0.0; 3],
50];
51
52pub const CQ: [[f64; 3]; 4] = [
53    [-0.500043713, 2.000209381, 3.135827145],
54    [6.531869153, -6.78386584, 7.247588801],
55    [-16.01477983, 20.38324603, 3.075947834],
56    [14.42597018, -10.89598394, 0.0],
57];
58
59// Dipole-Quadrupole parameters
60pub const ADQ: [[f64; 3]; 4] = [
61    [0.697094963, -0.673459279, 0.670340770],
62    [-0.633554144, -1.425899106, -4.338471826],
63    [2.945509028, 4.19441392, 7.234168360],
64    [-1.467027314, 1.0266216, 0.0],
65];
66
67pub const BDQ: [[f64; 3]; 4] = [
68    [-0.484038322, 0.67651011, -1.167560146],
69    [1.970405465, -3.013867512, 2.13488432],
70    [-2.118572671, 0.46742656, 0.0],
71    [0.0; 3],
72];
73
74pub const CDQ: [[f64; 2]; 3] = [
75    [0.795009692, -2.099579397],
76    [3.386863396, -5.941376392],
77    [0.475106328, -0.178820384],
78];
79
80pub const PI_SQ_43: f64 = 4.0 * PI * FRAC_PI_3;
81
82pub struct MeanSegmentNumbers {
83    pub mij1: Array2<f64>,
84    pub mij2: Array2<f64>,
85    pub mijk1: Array3<f64>,
86    pub mijk2: Array3<f64>,
87}
88
89pub enum Multipole {
90    Dipole,
91    Quadrupole,
92}
93
94impl MeanSegmentNumbers {
95    pub fn new(parameters: &PcSaftParameters, polarity: Multipole) -> Self {
96        let (npoles, comp) = match polarity {
97            Multipole::Dipole => (parameters.ndipole, &parameters.dipole_comp),
98            Multipole::Quadrupole => (parameters.nquadpole, &parameters.quadpole_comp),
99        };
100
101        let mut mi;
102        let mut mj;
103        let mut mk;
104        let mut mij;
105        let mut mijk;
106        let mut mij1 = Array2::zeros((npoles, npoles));
107        let mut mij2 = Array2::zeros((npoles, npoles));
108        let mut mijk1 = Array3::zeros((npoles, npoles, npoles));
109        let mut mijk2 = Array3::zeros((npoles, npoles, npoles));
110        for i in 0..npoles {
111            let dci = comp[i];
112            mi = parameters.m[dci].min(2.0);
113            mij1[[i, i]] = (mi - 1.0) / mi;
114            mij2[[i, i]] = mij1[[i, i]] * (mi - 2.0) / mi;
115
116            mijk1[[i, i, i]] = mij1[[i, i]];
117            mijk2[[i, i, i]] = mij2[[i, i]];
118            for j in i + 1..npoles {
119                let dcj = comp[j];
120                mj = parameters.m[dcj].min(2.0);
121                mij = (mi * mj).sqrt();
122                mij1[[i, j]] = (mij - 1.0) / mij;
123                mij2[[i, j]] = mij1[[i, j]] * (mij - 2.0) / mij;
124                mijk = (mi * mi * mj).cbrt();
125                mijk1[[i, i, j]] = (mijk - 1.0) / mijk;
126                mijk2[[i, i, j]] = mijk1[[i, i, j]] * (mijk - 2.0) / mijk;
127                mijk = (mi * mj * mj).cbrt();
128                mijk1[[i, j, j]] = (mijk - 1.0) / mijk;
129                mijk2[[i, j, j]] = mijk1[[i, j, j]] * (mijk - 2.0) / mijk;
130                for k in j + 1..npoles {
131                    let dck = comp[k];
132                    mk = parameters.m[dck].min(2.0);
133                    mijk = (mi * mj * mk).cbrt();
134                    mijk1[[i, j, k]] = (mijk - 1.0) / mijk;
135                    mijk2[[i, j, k]] = mijk1[[i, j, k]] * (mijk - 2.0) / mijk;
136                }
137            }
138        }
139        Self {
140            mij1,
141            mij2,
142            mijk1,
143            mijk2,
144        }
145    }
146}
147
148fn pair_integral_ij<D: DualNum<f64>>(
149    mij1: f64,
150    mij2: f64,
151    eta: D,
152    a: &[[f64; 3]],
153    b: &[[f64; 3]],
154    eps_ij_t: D,
155) -> D {
156    let eta2 = eta * eta;
157    let etas = [D::one(), eta, eta2, eta2 * eta, eta2 * eta2];
158    (0..a.len())
159        .map(|i| {
160            etas[i]
161                * (eps_ij_t * (b[i][0] + mij1 * b[i][1] + mij2 * b[i][2])
162                    + a[i][0]
163                    + mij1 * a[i][1]
164                    + mij2 * a[i][2])
165        })
166        .sum()
167}
168
169fn triplet_integral_ijk<D: DualNum<f64>>(mijk1: f64, mijk2: f64, eta: D, c: &[[f64; 3]]) -> D {
170    let eta2 = eta * eta;
171    let etas = [D::one(), eta, eta2, eta2 * eta];
172    (0..c.len())
173        .map(|i| etas[i] * (c[i][0] + mijk1 * c[i][1] + mijk2 * c[i][2]))
174        .sum()
175}
176
177fn triplet_integral_ijk_dq<D: DualNum<f64>>(mijk: f64, eta: D, c: &[[f64; 2]]) -> D {
178    let etas = [D::one(), eta, eta * eta];
179    (0..c.len())
180        .map(|i| etas[i] * (c[i][0] + mijk * c[i][1]))
181        .sum()
182}
183
184pub struct Dipole {
185    pub parameters: Rc<PcSaftParameters>,
186}
187
188impl<D: DualNum<f64>> HelmholtzEnergyDual<D> for Dipole {
189    fn helmholtz_energy(&self, state: &StateHD<D>) -> D {
190        let m = MeanSegmentNumbers::new(&self.parameters, Multipole::Dipole);
191        let p = &self.parameters;
192
193        let t_inv = state.temperature.inv();
194        let eps_ij_t = p.e_k_ij.mapv(|v| t_inv * v);
195        let sig_ij_3 = p.sigma_ij.mapv(|v| v.powi(3));
196        let mu2_term: Array1<D> = p
197            .dipole_comp
198            .iter()
199            .map(|&i| t_inv * sig_ij_3[[i, i]] * p.epsilon_k[i] * p.mu2[i])
200            .collect();
201
202        let rho = &state.partial_density;
203        let r = p.hs_diameter(state.temperature) * 0.5;
204        let eta = (rho * &p.m * &r * &r * &r).sum() * 4.0 * FRAC_PI_3;
205
206        let mut phi2 = D::zero();
207        let mut phi3 = D::zero();
208        for i in 0..p.ndipole {
209            let di = p.dipole_comp[i];
210            phi2 -= rho[di]
211                * rho[di]
212                * mu2_term[i]
213                * mu2_term[i]
214                * pair_integral_ij(
215                    m.mij1[[i, i]],
216                    m.mij2[[i, i]],
217                    eta,
218                    &AD,
219                    &BD,
220                    eps_ij_t[[di, di]],
221                )
222                / sig_ij_3[[di, di]];
223            phi3 -= rho[di]
224                * rho[di]
225                * rho[di]
226                * mu2_term[i]
227                * mu2_term[i]
228                * mu2_term[i]
229                * triplet_integral_ijk(m.mijk1[[i, i, i]], m.mijk2[[i, i, i]], eta, &CD)
230                / sig_ij_3[[di, di]];
231            for j in (i + 1)..p.ndipole {
232                let dj = p.dipole_comp[j];
233                phi2 -= rho[di]
234                    * rho[dj]
235                    * mu2_term[i]
236                    * mu2_term[j]
237                    * pair_integral_ij(
238                        m.mij1[[i, j]],
239                        m.mij2[[i, j]],
240                        eta,
241                        &AD,
242                        &BD,
243                        eps_ij_t[[di, dj]],
244                    )
245                    / sig_ij_3[[di, dj]]
246                    * 2.0;
247                phi3 -= rho[di] * rho[di] * rho[dj] * mu2_term[i] * mu2_term[i] * mu2_term[j]
248                    / (p.sigma_ij[[di, di]] * p.sigma_ij[[di, dj]] * p.sigma_ij[[di, dj]])
249                    * triplet_integral_ijk(m.mijk1[[i, i, j]], m.mijk2[[i, i, j]], eta, &CD)
250                    * 3.0;
251                phi3 -= rho[di] * rho[dj] * rho[dj] * mu2_term[i] * mu2_term[j] * mu2_term[j]
252                    / (p.sigma_ij[[di, dj]] * p.sigma_ij[[di, dj]] * p.sigma_ij[[dj, dj]])
253                    * triplet_integral_ijk(m.mijk1[[i, j, j]], m.mijk2[[i, j, j]], eta, &CD)
254                    * 3.0;
255                for k in (j + 1)..p.ndipole {
256                    let dk = p.dipole_comp[k];
257                    phi3 -= rho[di] * rho[dj] * rho[dk] * mu2_term[i] * mu2_term[j] * mu2_term[k]
258                        / (p.sigma_ij[[di, dj]] * p.sigma_ij[[di, dk]] * p.sigma_ij[[dj, dk]])
259                        * triplet_integral_ijk(m.mijk1[[i, j, k]], m.mijk2[[i, j, k]], eta, &CD)
260                        * 6.0;
261                }
262            }
263        }
264        phi2 *= PI;
265        phi3 *= PI_SQ_43;
266        let mut result = phi2 * phi2 / (phi2 - phi3) * state.volume;
267        if result.re().is_nan() {
268            result = phi2 * state.volume
269        }
270        result
271    }
272}
273
274impl fmt::Display for Dipole {
275    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
276        write!(f, "Dipole")
277    }
278}
279
280pub struct Quadrupole {
281    pub parameters: Rc<PcSaftParameters>,
282}
283
284impl<D: DualNum<f64>> HelmholtzEnergyDual<D> for Quadrupole {
285    fn helmholtz_energy(&self, state: &StateHD<D>) -> D {
286        let m = MeanSegmentNumbers::new(&self.parameters, Multipole::Quadrupole);
287        let p = &self.parameters;
288
289        let t_inv = state.temperature.inv();
290        let eps_ij_t = p.e_k_ij.mapv(|v| t_inv * v);
291        let sig_ij_3 = p.sigma_ij.mapv(|v| v.powi(3));
292        let q2_term: Array1<D> = p
293            .quadpole_comp
294            .iter()
295            .map(|&i| t_inv * p.sigma[i].powi(5) * p.epsilon_k[i] * p.q2[i])
296            .collect();
297
298        let rho = &state.partial_density;
299        let r = p.hs_diameter(state.temperature) * 0.5;
300        let eta = (rho * &p.m * &r * &r * &r).sum() * 4.0 * FRAC_PI_3;
301
302        let mut phi2 = D::zero();
303        let mut phi3 = D::zero();
304        for i in 0..p.nquadpole {
305            let di = p.quadpole_comp[i];
306            phi2 -= (rho[di]
307                * rho[di]
308                * q2_term[i]
309                * q2_term[i]
310                * pair_integral_ij(
311                    m.mij1[[i, i]],
312                    m.mij2[[i, i]],
313                    eta,
314                    &AQ,
315                    &BQ,
316                    eps_ij_t[[di, di]],
317                ))
318                / p.sigma_ij[[di, di]].powi(7);
319            phi3 += (rho[di]
320                * rho[di]
321                * rho[di]
322                * q2_term[i]
323                * q2_term[i]
324                * q2_term[i]
325                * triplet_integral_ijk(m.mijk1[[i, i, i]], m.mijk2[[i, i, i]], eta, &CQ))
326                / sig_ij_3[[di, di]].powi(3);
327            for j in (i + 1)..p.nquadpole {
328                let dj = p.quadpole_comp[j];
329                phi2 -= (rho[di]
330                    * rho[dj]
331                    * q2_term[i]
332                    * q2_term[j]
333                    * pair_integral_ij(
334                        m.mij1[[i, j]],
335                        m.mij2[[i, j]],
336                        eta,
337                        &AQ,
338                        &BQ,
339                        eps_ij_t[[di, dj]],
340                    ))
341                    / p.sigma_ij[[di, di]].powi(7)
342                    * 2.0;
343                phi3 += rho[di] * rho[di] * rho[dj] * q2_term[i] * q2_term[i] * q2_term[j]
344                    / (sig_ij_3[[di, di]] * sig_ij_3[[di, dj]] * sig_ij_3[[di, dj]])
345                    * triplet_integral_ijk(m.mijk1[[i, i, j]], m.mijk2[[i, i, j]], eta, &CQ)
346                    * 3.0;
347                phi3 += rho[di] * rho[dj] * rho[dj] * q2_term[i] * q2_term[j] * q2_term[j]
348                    / (sig_ij_3[[di, dj]] * sig_ij_3[[di, dj]] * sig_ij_3[[dj, dj]])
349                    * triplet_integral_ijk(m.mijk1[[i, j, j]], m.mijk2[[i, j, j]], eta, &CQ)
350                    * 3.0;
351                for k in (j + 1)..p.nquadpole {
352                    let dk = p.quadpole_comp[k];
353                    phi3 += rho[di] * rho[dj] * rho[dk] * q2_term[i] * q2_term[j] * q2_term[k]
354                        / (sig_ij_3[[di, dj]] * sig_ij_3[[di, dk]] * sig_ij_3[[dj, dk]])
355                        * triplet_integral_ijk(m.mijk1[[i, j, k]], m.mijk2[[i, j, k]], eta, &CQ)
356                        * 6.0;
357                }
358            }
359        }
360
361        phi2 *= PI * 0.5625;
362        phi3 *= PI * PI * 0.5625;
363        let mut result = phi2 * phi2 / (phi2 - phi3) * state.volume;
364        if result.re().is_nan() {
365            result = phi2 * state.volume
366        }
367        result
368    }
369}
370
371impl fmt::Display for Quadrupole {
372    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
373        write!(f, "Quadrupole")
374    }
375}
376
377#[derive(Clone, Copy)]
378pub enum DQVariants {
379    DQ35,
380    DQ44,
381}
382
383pub struct DipoleQuadrupole {
384    pub parameters: Rc<PcSaftParameters>,
385    pub variant: DQVariants,
386}
387
388impl<D: DualNum<f64>> HelmholtzEnergyDual<D> for DipoleQuadrupole {
389    fn helmholtz_energy(&self, state: &StateHD<D>) -> D {
390        let p = &self.parameters;
391
392        let t_inv = state.temperature.inv();
393        let eps_ij_t = p.e_k_ij.mapv(|v| t_inv * v);
394
395        let q2_term: Array1<D> = p
396            .quadpole_comp
397            .iter()
398            .map(|&i| t_inv * p.sigma[i].powi(4) * p.epsilon_k[i] * p.q2[i])
399            .collect();
400        let mu2_term: Array1<D> = p
401            .dipole_comp
402            .iter()
403            .map(|&i| t_inv * p.sigma[i].powi(4) * p.epsilon_k[i] * p.mu2[i])
404            .collect();
405
406        let rho = &state.partial_density;
407        let r = p.hs_diameter(state.temperature) * 0.5;
408        let eta = (rho * &p.m * &r * &r * &r).sum() * 4.0 * FRAC_PI_3;
409
410        // mean segment number
411        let mut mdq1 = Array2::zeros((p.ndipole, p.nquadpole));
412        let mut mdq2 = Array2::zeros((p.ndipole, p.nquadpole));
413        let mut mdqd = Array3::zeros((p.ndipole, p.nquadpole, p.ndipole));
414        let mut mdqq = Array3::zeros((p.ndipole, p.nquadpole, p.nquadpole));
415        for i in 0..p.ndipole {
416            let di = p.dipole_comp[i];
417            let mi = p.m[di].min(2.0);
418            for j in 0..p.nquadpole {
419                let qj = p.quadpole_comp[j];
420                let mj = p.m[qj].min(2.0);
421                let m = (mi * mj).sqrt();
422                mdq1[[i, j]] = (m - 1.0) / m;
423                mdq2[[i, j]] = mdq1[[i, j]] * (m - 2.0) / m;
424                for k in 0..p.ndipole {
425                    let dk = p.dipole_comp[k];
426                    let mk = p.m[dk].min(2.0);
427                    let m = (mi * mj * mk).cbrt();
428                    mdqd[[i, j, k]] = (m - 1.0) / m;
429                }
430                for k in 0..p.nquadpole {
431                    let qk = p.quadpole_comp[k];
432                    let mk = p.m[qk].min(2.0);
433                    let m = (mi * mj * mk).cbrt();
434                    mdqq[[i, j, k]] = (m - 1.0) / m;
435                }
436            }
437        }
438
439        let mut phi2 = D::zero();
440        let mut phi3 = D::zero();
441        for i in 0..p.ndipole {
442            for j in 0..p.nquadpole {
443                let di = p.dipole_comp[i];
444                let qj = p.quadpole_comp[j];
445                let mu2_q2_term = match self.variant {
446                    DQVariants::DQ35 => mu2_term[i] / p.sigma[di] * q2_term[j] * p.sigma[qj],
447                    DQVariants::DQ44 => mu2_term[i] * q2_term[j],
448                };
449                phi2 -= rho[di] * rho[qj] * mu2_q2_term / p.sigma_ij[[di, qj]].powi(5)
450                    * pair_integral_ij(
451                        mdq1[[i, j]],
452                        mdq2[[i, j]],
453                        eta,
454                        &ADQ,
455                        &BDQ,
456                        eps_ij_t[[di, qj]],
457                    );
458                for k in 0..p.ndipole {
459                    let dk = p.dipole_comp[k];
460                    phi3 += rho[di] * rho[qj] * rho[dk] * mu2_term[i] * q2_term[j] * mu2_term[k]
461                        / (p.sigma_ij[[di, qj]] * p.sigma_ij[[di, dk]] * p.sigma_ij[[qj, dk]])
462                            .powi(2)
463                        * triplet_integral_ijk_dq(mdqd[[i, j, k]], eta, &CDQ);
464                }
465                for k in 0..p.nquadpole {
466                    let qk = p.quadpole_comp[k];
467                    phi3 += rho[di] * rho[qj] * rho[qk] * mu2_term[i] * q2_term[j] * q2_term[k]
468                        / (p.sigma_ij[[di, qj]] * p.sigma_ij[[di, qk]] * p.sigma_ij[[qj, qk]])
469                            .powi(2)
470                        * ALPHA
471                        * triplet_integral_ijk_dq(mdqq[[i, j, k]], eta, &CDQ);
472                }
473            }
474        }
475
476        phi2 *= PI * 2.25;
477        phi3 *= PI * PI;
478        let mut result = phi2 * phi2 / (phi2 - phi3) * state.volume;
479        if result.re().is_nan() {
480            result = phi2 * state.volume
481        }
482        result
483    }
484}
485
486impl fmt::Display for DipoleQuadrupole {
487    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
488        write!(f, "DipoleQuadrupole")
489    }
490}
491
492#[cfg(test)]
493mod tests {
494    use super::*;
495    use crate::eos::dispersion::Dispersion;
496    use crate::parameters::utils::{carbon_dioxide_parameters, dme_co2_parameters, dme_parameters};
497    use approx::assert_relative_eq;
498    use feos_core::StateHD;
499
500    #[test]
501    fn test_dipolar_contribution() {
502        let dp = Dipole {
503            parameters: Rc::new(dme_parameters()),
504        };
505        let t = 350.0;
506        let v = 1000.0;
507        let n = 1.0;
508        let s = StateHD::new(t, v, arr1(&[n]));
509        let a = dp.helmholtz_energy(&s);
510        assert_relative_eq!(a, -1.40501033595417E-002, epsilon = 1e-6);
511    }
512
513    #[test]
514    fn test_quadrupolar_contribution() {
515        let qp = Quadrupole {
516            parameters: Rc::new(carbon_dioxide_parameters()),
517        };
518        let t = 350.0;
519        let v = 1000.0;
520        let n = 1.0;
521        let s = StateHD::new(t, v, arr1(&[n]));
522        let a = qp.helmholtz_energy(&s);
523        assert_relative_eq!(a, -4.38559558854186E-002, epsilon = 1e-6);
524    }
525
526    #[test]
527    fn test_dipolar_quadrupolar_contribution() {
528        let dp = Dipole {
529            parameters: Rc::new(dme_co2_parameters()),
530        };
531        let qp = Quadrupole {
532            parameters: Rc::new(dme_co2_parameters()),
533        };
534        // let dpqp = DipoleQuadrupole {
535        //     parameters: Rc::new(dme_co2_parameters()),
536        //     variant: DQVariants::DQ35,
537        // };
538        let disp = Dispersion {
539            parameters: Rc::new(dme_co2_parameters()),
540        };
541        let t = 350.0;
542        let v = 1000.0;
543        let n_dme = 1.0;
544        let n_co2 = 1.0;
545        let s = StateHD::new(t, v, arr1(&[n_dme, n_co2]));
546        // let a_dpqp = dpqp.helmholtz_energy(&s);
547        let a_disp = disp.helmholtz_energy(&s);
548        let a_qp = qp.helmholtz_energy(&s);
549        let a_dp = dp.helmholtz_energy(&s);
550        assert_relative_eq!(a_disp, -1.6283622072860044, epsilon = 1e-6);
551        assert_relative_eq!(a_dp, -1.35361827881345E-002, epsilon = 1e-6);
552        assert_relative_eq!(a_qp, -4.20168059082731E-002, epsilon = 1e-6);
553        // assert_relative_eq!(a_dpqp, -2.2316252638709004E-002, epsilon = 1e-6);
554    }
555}
feos_pcsaft/eos/polar.rs

feos_pcsaft/eos/
polar.rs
