Repository navigation
Expand file tree
/
Copy pathdispersion.rs
More file actions
136 lines (122 loc) · 4.73 KB
/
Copy pathdispersion.rs
File metadata and controls
136 lines (122 loc) · 4.73 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
use super::polar::calculate_helmholtz_energy_density_polar;
use super::PcSaftParameters;
use crate::hard_sphere::HardSphereProperties;
use crate::pcsaft::eos::dispersion::{A0, A1, A2, B0, B1, B2};
use feos_core::EosError;
use feos_dft::{
FunctionalContributionDual, WeightFunction, WeightFunctionInfo, WeightFunctionShape,
};
use ndarray::*;
use num_dual::DualNum;
use std::f64::consts::{FRAC_PI_3, PI};
use std::fmt;
use std::sync::Arc;
/// psi Parameter for DFT (Sauer2017)
const PSI_DFT: f64 = 1.3862;
/// psi Parameter for pDGT (Rehner2018)
const PSI_PDGT: f64 = 1.3286;
#[derive(Clone)]
pub struct AttractiveFunctional {
parameters: Arc<PcSaftParameters>,
}
impl AttractiveFunctional {
pub fn new(parameters: Arc<PcSaftParameters>) -> Self {
Self { parameters }
}
}
fn att_weight_functions<N: DualNum<f64> + Copy + ScalarOperand>(
p: &PcSaftParameters,
psi: f64,
temperature: N,
) -> WeightFunctionInfo<N> {
let d = p.hs_diameter(temperature);
WeightFunctionInfo::new(Array1::from_shape_fn(d.len(), |i| i), false).add(
WeightFunction::new_scaled(d * psi, WeightFunctionShape::Theta),
false,
)
}
impl<N: DualNum<f64> + Copy + ScalarOperand> FunctionalContributionDual<N>
for AttractiveFunctional
{
fn weight_functions(&self, temperature: N) -> WeightFunctionInfo<N> {
att_weight_functions(&self.parameters, PSI_DFT, temperature)
}
fn weight_functions_pdgt(&self, temperature: N) -> WeightFunctionInfo<N> {
att_weight_functions(&self.parameters, PSI_PDGT, temperature)
}
fn calculate_helmholtz_energy_density(
&self,
temperature: N,
density: ArrayView2<N>,
) -> Result<Array1<N>, EosError> {
// auxiliary variables
let p = &self.parameters;
let n = p.m.len();
// temperature dependent segment radius
let r = p.hs_diameter(temperature) * 0.5;
// packing fraction
let eta = density
.outer_iter()
.zip((&r * &r * &r * &p.m * 4.0 * FRAC_PI_3).into_iter())
.fold(
Array::zeros(density.raw_dim().remove_axis(Axis(0))),
|acc: Array1<N>, (rho, r3m)| acc + &rho * r3m,
);
// mean segment number
let mut rhog = Array::zeros(eta.raw_dim());
let mut m_bar = Array::zeros(eta.raw_dim());
for (rhoi, &mi) in density.axis_iter(Axis(0)).zip(p.m.iter()) {
m_bar += &(&rhoi * mi);
rhog += &rhoi;
}
m_bar.iter_mut().zip(rhog.iter()).for_each(|(m, &r)| {
if r.re() > f64::EPSILON {
*m /= r
} else {
*m = N::one()
}
});
// mixture densities, crosswise interactions of all segments on all chains
let mut rho1mix: Array1<N> = Array::zeros(eta.raw_dim());
let mut rho2mix: Array1<N> = Array::zeros(eta.raw_dim());
for i in 0..n {
for j in 0..n {
let eps_ij_t = temperature.recip() * p.epsilon_k_ij[(i, j)];
let sigma_ij_3 = p.sigma_ij[(i, j)].powi(3);
rho1mix = rho1mix
+ (&density.index_axis(Axis(0), i) * &density.index_axis(Axis(0), j))
.mapv(|x| x * (eps_ij_t * sigma_ij_3 * p.m[i] * p.m[j]));
rho2mix = rho2mix
+ (&density.index_axis(Axis(0), i) * &density.index_axis(Axis(0), j))
.mapv(|x| x * (eps_ij_t * eps_ij_t * sigma_ij_3 * p.m[i] * p.m[j]));
}
}
// I1, I2 and C1
let mut i1: Array1<N> = Array::zeros(eta.raw_dim());
let mut i2: Array1<N> = Array::zeros(eta.raw_dim());
let mut eta_i: Array1<N> = Array::ones(eta.raw_dim());
let m1 = (m_bar.clone() - 1.0) / &m_bar;
let m2 = (m_bar.clone() - 2.0) / &m_bar * &m1;
for i in 0..=6 {
i1 = i1 + (&m2 * A2[i] + &m1 * A1[i] + A0[i]) * &eta_i;
i2 = i2 + (&m2 * B2[i] + &m1 * B1[i] + B0[i]) * &eta_i;
eta_i = &eta_i * η
}
let c1 = Zip::from(&eta).and(&m_bar).map_collect(|&eta, &m| {
(m * (eta * 8.0 - eta.powi(2) * 2.0) / (eta - 1.0).powi(4)
+ (eta * (eta * (eta * (eta * 2.0 - 12.0) + 27.0) - 20.0))
/ ((eta - 1.0) * (eta - 2.0)).powi(2)
* (m - 1.0)
+ 1.0)
.recip()
});
// Helmholtz energy density
let phi_polar = calculate_helmholtz_energy_density_polar(p, temperature, density)?;
Ok((-rho1mix * i1 * 2.0 - rho2mix * m_bar * c1 * i2) * PI + phi_polar)
}
}
impl fmt::Display for AttractiveFunctional {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
write!(f, "Attractive functional")
}
}