feos_core/
cubic.rs

1//! Implementation of the Peng-Robinson equation of state.
2//!
3//! This module acts as a reference on how a simple equation
4//! of state - with a single contribution to the Helmholtz energy - can be implemented.
5//! The implementation closely follows the form of the equations given in
6//! [this wikipedia article](https://en.wikipedia.org/wiki/Cubic_equations_of_state#Peng%E2%80%93Robinson_equation_of_state).
7use crate::parameter::{Identifier, Parameters, PureRecord};
8use crate::{FeosError, FeosResult, Molarweight, ResidualDyn, StateHD, Subset};
9use nalgebra::{DMatrix, DVector};
10use num_dual::DualNum;
11use quantity::MolarWeight;
12use serde::{Deserialize, Serialize};
13use std::f64::consts::SQRT_2;
14
15const KB_A3: f64 = 13806490.0;
16
17/// Peng-Robinson parameters for a single substance.
18#[derive(Serialize, Deserialize, Debug, Clone)]
19pub struct PengRobinsonRecord {
20    /// critical temperature in Kelvin
21    tc: f64,
22    /// critical pressure in Pascal
23    pc: f64,
24    /// acentric factor
25    acentric_factor: f64,
26}
27
28impl PengRobinsonRecord {
29    /// Create a new pure substance record for the Peng-Robinson equation of state.
30    pub fn new(tc: f64, pc: f64, acentric_factor: f64) -> Self {
31        Self {
32            tc,
33            pc,
34            acentric_factor,
35        }
36    }
37}
38
39/// Peng-Robinson parameters for one ore more substances.
40pub type PengRobinsonParameters = Parameters<PengRobinsonRecord, f64, ()>;
41
42impl PengRobinsonParameters {
43    /// Build a simple parameter set without binary interaction parameters.
44    pub fn new_simple(
45        tc: &[f64],
46        pc: &[f64],
47        acentric_factor: &[f64],
48        molarweight: &[f64],
49    ) -> FeosResult<Self> {
50        if [pc.len(), acentric_factor.len(), molarweight.len()]
51            .iter()
52            .any(|&l| l != tc.len())
53        {
54            return Err(FeosError::IncompatibleParameters(String::from(
55                "each component has to have parameters.",
56            )));
57        }
58        let records = (0..tc.len())
59            .map(|i| {
60                let record = PengRobinsonRecord {
61                    tc: tc[i],
62                    pc: pc[i],
63                    acentric_factor: acentric_factor[i],
64                };
65                let id = Identifier::default();
66                PureRecord::new(id, molarweight[i], record)
67            })
68            .collect();
69        PengRobinsonParameters::new(records, vec![])
70    }
71}
72
73/// A simple version of the Peng-Robinson equation of state.
74pub struct PengRobinson {
75    /// Parameters
76    pub parameters: PengRobinsonParameters,
77    /// Critical temperature in Kelvin
78    tc: DVector<f64>,
79    a: DVector<f64>,
80    b: DVector<f64>,
81    /// Binary interaction parameter
82    k_ij: DMatrix<f64>,
83    kappa: DVector<f64>,
84}
85
86impl PengRobinson {
87    /// Create a new equation of state from a set of parameters.
88    pub fn new(parameters: PengRobinsonParameters) -> Self {
89        let [tc, pc, ac] = parameters.collate(|r| [r.tc, r.pc, r.acentric_factor]);
90        let [k_ij] = parameters.collate_binary(|&br| [br]);
91
92        let a = 0.45724 * KB_A3 * tc.component_mul(&tc).component_div(&pc);
93        let b = 0.07780 * KB_A3 * &tc.component_div(&pc);
94        let kappa = ac.map(|ac| 0.37464 + (1.54226 - 0.26992 * &ac) * ac);
95        Self {
96            parameters,
97            tc,
98            a,
99            b,
100            k_ij,
101            kappa,
102        }
103    }
104}
105
106impl ResidualDyn for PengRobinson {
107    fn components(&self) -> usize {
108        self.tc.len()
109    }
110
111    fn compute_max_density<D: DualNum<f64> + Copy>(&self, molefracs: &DVector<D>) -> D {
112        D::from(0.9) / molefracs.dot(&self.b.map(D::from))
113    }
114
115    fn reduced_helmholtz_energy_density_contributions<D: DualNum<f64> + Copy>(
116        &self,
117        state: &StateHD<D>,
118    ) -> Vec<(&'static str, D)> {
119        let density = state.partial_density.sum();
120        let x = &state.molefracs;
121        let ak = &self
122            .tc
123            .map(|tc| D::one() - (state.temperature / tc).sqrt())
124            .component_mul(&self.kappa.map(D::from))
125            .map(|x| (x + 1.0).powi(2))
126            .component_mul(&self.a.map(D::from));
127
128        // Mixing rules
129        let mut ak_mix = D::zero();
130        for i in 0..ak.len() {
131            for j in 0..ak.len() {
132                ak_mix += (ak[i] * ak[j]).sqrt() * (x[i] * x[j] * (1.0 - self.k_ij[(i, j)]));
133            }
134        }
135        let b = x.dot(&self.b.map(D::from));
136
137        // Helmholtz energy
138        let v = density.recip();
139        let f = density
140            * ((v / (v - b)).ln()
141                - ak_mix / (b * SQRT_2 * 2.0 * state.temperature)
142                    * ((v + b * (1.0 + SQRT_2)) / (v + b * (1.0 - SQRT_2))).ln());
143        vec![("Peng Robinson", f)]
144    }
145}
146
147impl Subset for PengRobinson {
148    fn subset(&self, component_list: &[usize]) -> Self {
149        Self::new(self.parameters.subset(component_list))
150    }
151}
152
153impl Molarweight for PengRobinson {
154    fn molar_weight(&self) -> MolarWeight<DVector<f64>> {
155        self.parameters.molar_weight.clone()
156    }
157}
158
159#[cfg(test)]
160mod tests {
161    use super::*;
162    use crate::parameter::PureRecord;
163    use crate::state::{Contributions, State};
164    use crate::{FeosResult, SolverOptions, Verbosity};
165    use approx::*;
166    use quantity::{KELVIN, PASCAL};
167
168    fn pure_record_vec() -> Vec<PureRecord<PengRobinsonRecord, ()>> {
169        let records = r#"[
170            {
171                "identifier": {
172                    "cas": "74-98-6",
173                    "name": "propane",
174                    "iupac_name": "propane",
175                    "smiles": "CCC",
176                    "inchi": "InChI=1/C3H8/c1-3-2/h3H2,1-2H3",
177                    "formula": "C3H8"
178                },
179                "tc": 369.96,
180                "pc": 4250000.0,
181                "acentric_factor": 0.153,
182                "molarweight": 44.0962
183            },
184            {
185                "identifier": {
186                    "cas": "106-97-8",
187                    "name": "butane",
188                    "iupac_name": "butane",
189                    "smiles": "CCCC",
190                    "inchi": "InChI=1/C4H10/c1-3-4-2/h3-4H2,1-2H3",
191                    "formula": "C4H10"
192                },
193                "tc": 425.2,
194                "pc": 3800000.0,
195                "acentric_factor": 0.199,
196                "molarweight": 58.123
197            }
198        ]"#;
199        serde_json::from_str(records).expect("Unable to parse json.")
200    }
201
202    #[test]
203    fn peng_robinson() -> FeosResult<()> {
204        let mixture = pure_record_vec();
205        let propane = mixture[0].clone();
206        let tc = propane.model_record.tc;
207        let pc = propane.model_record.pc;
208        let parameters = PengRobinsonParameters::new_pure(propane)?;
209        let pr = PengRobinson::new(parameters);
210        let options = SolverOptions::new().verbosity(Verbosity::Iter);
211        let cp = State::critical_point(&&pr, None, None, None, options)?;
212        println!("{} {}", cp.temperature, cp.pressure(Contributions::Total));
213        assert_relative_eq!(cp.temperature, tc * KELVIN, max_relative = 1e-4);
214        assert_relative_eq!(
215            cp.pressure(Contributions::Total),
216            pc * PASCAL,
217            max_relative = 1e-4
218        );
219        Ok(())
220    }
221}
feos_core/cubic.rs

feos_core/
cubic.rs
