feos_ad/eos/
ideal_gas.rs

1use crate::{IdealGasAD, ParametersAD};
2use num_dual::DualNum;
3use std::collections::HashMap;
4
5const RGAS: f64 = 6.022140857 * 1.38064852;
6const T0: f64 = 298.15;
7const T0_2: f64 = 298.15 * 298.15;
8const T0_3: f64 = T0 * T0_2;
9const T0_4: f64 = T0_2 * T0_2;
10const T0_5: f64 = T0 * T0_4;
11const P0: f64 = 1.0e5;
12const A3: f64 = 1e-30;
13const KB: f64 = 1.38064852e-23;
14
15const A: [f64; 22] = [
16    19.5, -0.909, -23.0, -66.2, 23.6, -8.0, -28.0, 32.37, -6.03, -20.5, -6.03, -20.5, -2.14, -8.25,
17    30.9, 6.45, 45.0, 24.591, 24.1, 24.5, 25.7, 26.9,
18];
19const B: [f64; 22] = [
20    -0.00808, 0.095, 0.204, 0.427, -0.0381, 0.105, 0.208, -0.007, 0.0854, 0.162, 0.0854, 0.162,
21    0.0574, 0.101, -0.0336, 0.067, -0.07128, 0.0318, 0.0427, 0.0402, -0.0691, -0.0412,
22];
23const C: [f64; 22] = [
24    0.000153, -5.44e-05, -0.000265, -0.000641, 0.000172, -9.63e-05, -0.000306, 0.00010267, -8e-06,
25    -0.00016, -8e-06, -0.00016, -1.64e-06, -0.000142, 0.00016, -3.57e-05, 0.000264, 5.66e-05,
26    8.04e-05, 4.02e-05, 0.000177, 0.000164,
27];
28const D: [f64; 22] = [
29    -9.67e-08, 1.19e-08, 1.2e-07, 3.01e-07, -1.03e-07, 3.56e-08, 1.46e-07, -6.641e-08, -1.8e-08,
30    6.24e-08, -1.8e-08, 6.24e-08, -1.59e-08, 6.78e-08, -9.88e-08, 2.86e-09, -1.515e-07, -4.29e-08,
31    -6.87e-08, -4.52e-08, -9.88e-08, -9.76e-08,
32];
33const GROUPS: [&str; 22] = [
34    "CH3", "CH2", ">CH", ">C<", "=CH2", "=CH", "=C<", "C≡CH", "CH2_hex", "CH_hex", "CH2_pent",
35    "CH_pent", "CH_arom", "C_arom", "CH=O", ">C=O", "OCH3", "OCH2", "HCOO", "COO", "OH", "NH2",
36];
37
38/// The GC method for the ideal gas heat capacity by Joback & Reid.
39#[derive(Clone, Copy)]
40pub struct Joback(pub [f64; 5]);
41
42impl Joback {
43    pub fn from_groups<D: DualNum<f64> + Copy>(group_counts: [D; 22]) -> [D; 5] {
44        let a: D = A.into_iter().zip(group_counts).map(|(a, g)| g * a).sum();
45        let b: D = B.into_iter().zip(group_counts).map(|(b, g)| g * b).sum();
46        let c: D = C.into_iter().zip(group_counts).map(|(c, g)| g * c).sum();
47        let d: D = D.into_iter().zip(group_counts).map(|(d, g)| g * d).sum();
48
49        [a - 37.93, b + 0.21, c - 3.91e-4, d + 2.06e-7, D::zero()]
50    }
51
52    pub fn from_group_counts<D: DualNum<f64> + Copy>(group_counts: &HashMap<&str, D>) -> [D; 5] {
53        Self::from_groups(GROUPS.map(|g| *group_counts.get(g).unwrap_or(&D::zero())))
54    }
55}
56
57impl ParametersAD for Joback {
58    type Parameters<D: DualNum<f64> + Copy> = [D; 5];
59
60    fn params<D: DualNum<f64> + Copy>(&self) -> [D; 5] {
61        self.0.map(D::from)
62    }
63
64    fn params_from_inner<D: DualNum<f64> + Copy, D2: DualNum<f64, Inner = D> + Copy>(
65        parameters: &[D; 5],
66    ) -> [D2; 5] {
67        parameters.map(D2::from_inner)
68    }
69}
70
71impl IdealGasAD for Joback {
72    const IDEAL_GAS: &str = "Joback";
73
74    fn ln_lambda3_dual<D: DualNum<f64> + Copy>(
75        parameters: &Self::Parameters<D>,
76        temperature: D,
77    ) -> D {
78        let [a, b, c, d, e] = *parameters;
79        let t = temperature;
80        let t2 = t * t;
81        let t3 = t2 * t;
82        let t4 = t2 * t2;
83        let f = (temperature * KB / (P0 * A3)).ln();
84        let h = (t2 - T0_2) * 0.5 * b
85            + (t3 - T0_3) * c / 3.0
86            + (t4 - T0_4) * 0.25 * d
87            + (t4 * t - T0_5) * 0.2 * e
88            + (t - T0) * a;
89        let s = (t - T0) * b
90            + (t2 - T0_2) * 0.5 * c
91            + (t3 - T0_3) * d / 3.0
92            + (t4 - T0_4) * 0.25 * e
93            + (t / T0).ln() * a;
94        (h - t * s) / (t * RGAS) + f
95    }
96}
97
98#[cfg(test)]
99pub mod test {
100    use super::Joback as JobackAD;
101    use crate::{EquationOfStateAD, ParametersAD, ResidualHelmholtzEnergy, TotalHelmholtzEnergy};
102    use approx::assert_relative_eq;
103    use feos::ideal_gas::{Joback, JobackRecord};
104    use feos_core::parameter::Parameter;
105    use feos_core::{Contributions::IdealGas, EosResult, EquationOfState, ReferenceSystem, State};
106    use nalgebra::SVector;
107    use ndarray::arr1;
108    use num_dual::DualNum;
109    use quantity::{KELVIN, KILO, METER, MOL};
110    use std::sync::Arc;
111
112    pub fn joback() -> EosResult<(JobackAD, Arc<Joback>)> {
113        let a = 1.5;
114        let b = 3.4e-2;
115        let c = 180.0e-4;
116        let d = 2.2e-6;
117        let e = 0.03e-8;
118        let eos = Arc::new(Joback::from_model_records(vec![JobackRecord::new(
119            a, b, c, d, e,
120        )])?);
121        let params = [a, b, c, d, e];
122        let eos_ad = JobackAD(params);
123        Ok((eos_ad, eos))
124    }
125
126    struct NoResidual;
127    impl ParametersAD for NoResidual {
128        type Parameters<D: DualNum<f64> + Copy> = [D; 0];
129
130        fn params<D: DualNum<f64> + Copy>(&self) -> Self::Parameters<D> {
131            []
132        }
133
134        fn params_from_inner<D: DualNum<f64> + Copy, D2: DualNum<f64, Inner = D> + Copy>(
135            _: &Self::Parameters<D>,
136        ) -> Self::Parameters<D2> {
137            []
138        }
139    }
140    impl<const N: usize> ResidualHelmholtzEnergy<N> for NoResidual {
141        const RESIDUAL: &str = "No residual";
142
143        fn compute_max_density(&self, _: &SVector<f64, N>) -> f64 {
144            1.0
145        }
146
147        fn residual_helmholtz_energy_density<D: DualNum<f64> + Copy>(
148            _: &Self::Parameters<D>,
149            _: D,
150            _: &SVector<D, N>,
151        ) -> D {
152            D::zero()
153        }
154    }
155
156    type JobackEos = EquationOfStateAD<JobackAD, NoResidual, 1>;
157
158    #[test]
159    fn test_joback() -> EosResult<()> {
160        let (joback_ad, joback) = joback()?;
161        let eos = Arc::new(EquationOfState::ideal_gas(joback));
162        let eos_ad = EquationOfStateAD::new([joback_ad], NoResidual);
163
164        let temperature = 300.0 * KELVIN;
165        let volume = 2.3 * METER * METER * METER;
166        let moles = arr1(&[1.3]) * KILO * MOL;
167
168        let state = State::new_nvt(&eos, temperature, volume, &moles)?;
169        let a_feos = state.molar_helmholtz_energy(IdealGas);
170        let mu_feos = state.chemical_potential(IdealGas);
171        let p_feos = state.pressure(IdealGas);
172        let s_feos = state.molar_entropy(IdealGas);
173        let h_feos = state.molar_enthalpy(IdealGas);
174
175        let total_moles = moles.sum();
176        let t = temperature.to_reduced();
177        let v = (volume / total_moles).to_reduced();
178        let x = SVector::from_fn(|i, _| moles.get(i).convert_into(total_moles));
179        let a_ad = JobackEos::molar_helmholtz_energy(&eos_ad.params(), t, v, &x);
180        let mu_ad = JobackEos::chemical_potential(&eos_ad.params(), t, v, &x);
181        let p_ad = JobackEos::pressure(&eos_ad.params(), t, v, &x);
182        let s_ad = JobackEos::molar_entropy(&eos_ad.params(), t, v, &x);
183        let h_ad = JobackEos::molar_enthalpy(&eos_ad.params(), t, v, &x);
184
185        println!("\nMolar Helmholtz energy:\n{}", a_feos.to_reduced());
186        println!("{a_ad}");
187        assert_relative_eq!(a_feos.to_reduced(), a_ad, max_relative = 1e-14);
188
189        println!("\nChemical potential:\n{}", mu_feos.get(0).to_reduced());
190        println!("{}", mu_ad[0]);
191        assert_relative_eq!(mu_feos.get(0).to_reduced(), mu_ad[0], max_relative = 1e-14);
192
193        println!("\nPressure:\n{}", p_feos.to_reduced());
194        println!("{p_ad}");
195        assert_relative_eq!(p_feos.to_reduced(), p_ad, max_relative = 1e-14);
196
197        println!("\nMolar entropy:\n{}", s_feos.to_reduced());
198        println!("{s_ad}");
199        assert_relative_eq!(s_feos.to_reduced(), s_ad, max_relative = 1e-14);
200
201        println!("\nMolar enthalpy:\n{}", h_feos.to_reduced());
202        println!("{h_ad}");
203        assert_relative_eq!(h_feos.to_reduced(), h_ad, max_relative = 1e-14);
204
205        Ok(())
206    }
207}
feos_ad/eos/ideal_gas.rs

feos_ad/eos/
ideal_gas.rs
