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::equation_of_state::{Components, Molarweight, Residual};
8use crate::parameter::{Identifier, Parameter, ParameterError, PureRecord};
9use crate::state::StateHD;
10use ndarray::{Array1, Array2, ScalarOperand};
11use num_dual::DualNum;
12use quantity::{MolarWeight, GRAM, MOL};
13use serde::{Deserialize, Serialize};
14use std::f64::consts::SQRT_2;
15use std::fmt;
16use std::sync::Arc;
17
18const KB_A3: f64 = 13806490.0;
19
20/// Peng-Robinson parameters for a single substance.
21#[derive(Serialize, Deserialize, Debug, Clone, Default)]
22pub struct PengRobinsonRecord {
23    /// critical temperature in Kelvin
24    tc: f64,
25    /// critical pressure in Pascal
26    pc: f64,
27    /// acentric factor
28    acentric_factor: f64,
29}
30
31impl PengRobinsonRecord {
32    /// Create a new pure substance record for the Peng-Robinson equation of state.
33    pub fn new(tc: f64, pc: f64, acentric_factor: f64) -> Self {
34        Self {
35            tc,
36            pc,
37            acentric_factor,
38        }
39    }
40}
41
42impl std::fmt::Display for PengRobinsonRecord {
43    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
44        write!(f, "PengRobinsonRecord(tc={} K", self.tc)?;
45        write!(f, ", pc={} Pa", self.pc)?;
46        write!(f, ", acentric factor={}", self.acentric_factor)
47    }
48}
49
50/// Peng-Robinson parameters for one ore more substances.
51pub struct PengRobinsonParameters {
52    /// Critical temperature in Kelvin
53    tc: Array1<f64>,
54    a: Array1<f64>,
55    b: Array1<f64>,
56    /// Binary interaction parameter
57    k_ij: Array2<f64>,
58    kappa: Array1<f64>,
59    /// Molar weight in units of g/mol
60    molarweight: Array1<f64>,
61    /// List of pure component records
62    pure_records: Vec<PureRecord<PengRobinsonRecord>>,
63}
64
65impl std::fmt::Display for PengRobinsonParameters {
66    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
67        self.pure_records
68            .iter()
69            .try_for_each(|pr| writeln!(f, "{}", pr))?;
70        writeln!(f, "\nk_ij:\n{}", self.k_ij)
71    }
72}
73
74impl PengRobinsonParameters {
75    /// Build a simple parameter set without binary interaction parameters.
76    pub fn new_simple(
77        tc: &[f64],
78        pc: &[f64],
79        acentric_factor: &[f64],
80        molarweight: &[f64],
81    ) -> Result<Self, crate::parameter::ParameterError> {
82        if [pc.len(), acentric_factor.len(), molarweight.len()]
83            .iter()
84            .any(|&l| l != tc.len())
85        {
86            return Err(ParameterError::IncompatibleParameters(String::from(
87                "each component has to have parameters.",
88            )));
89        }
90        let records = (0..tc.len())
91            .map(|i| {
92                let record = PengRobinsonRecord {
93                    tc: tc[i],
94                    pc: pc[i],
95                    acentric_factor: acentric_factor[i],
96                };
97                let id = Identifier::default();
98                PureRecord::new(id, molarweight[i], record)
99            })
100            .collect();
101        PengRobinsonParameters::from_records(records, None)
102    }
103}
104
105impl Parameter for PengRobinsonParameters {
106    type Pure = PengRobinsonRecord;
107    type Binary = f64;
108
109    /// Creates parameters from pure component records.
110    fn from_records(
111        pure_records: Vec<PureRecord<Self::Pure>>,
112        binary_records: Option<Array2<Self::Binary>>,
113    ) -> Result<Self, ParameterError> {
114        let n = pure_records.len();
115
116        let mut tc = Array1::zeros(n);
117        let mut a = Array1::zeros(n);
118        let mut b = Array1::zeros(n);
119        let mut molarweight = Array1::zeros(n);
120        let mut kappa = Array1::zeros(n);
121
122        for (i, record) in pure_records.iter().enumerate() {
123            molarweight[i] = record.molarweight;
124            let r = &record.model_record;
125            tc[i] = r.tc;
126            a[i] = 0.45724 * r.tc.powi(2) * KB_A3 / r.pc;
127            b[i] = 0.07780 * r.tc * KB_A3 / r.pc;
128            kappa[i] = 0.37464 + (1.54226 - 0.26992 * r.acentric_factor) * r.acentric_factor;
129        }
130
131        let k_ij = binary_records.unwrap_or_else(|| Array2::zeros([n; 2]));
132
133        Ok(Self {
134            tc,
135            a,
136            b,
137            k_ij,
138            kappa,
139            molarweight,
140            pure_records,
141        })
142    }
143
144    fn records(&self) -> (&[PureRecord<PengRobinsonRecord>], Option<&Array2<f64>>) {
145        (&self.pure_records, Some(&self.k_ij))
146    }
147}
148
149/// A simple version of the Peng-Robinson equation of state.
150pub struct PengRobinson {
151    /// Parameters
152    parameters: Arc<PengRobinsonParameters>,
153}
154
155impl PengRobinson {
156    /// Create a new equation of state from a set of parameters.
157    pub fn new(parameters: Arc<PengRobinsonParameters>) -> Self {
158        Self { parameters }
159    }
160}
161
162impl fmt::Display for PengRobinson {
163    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
164        write!(f, "Peng Robinson")
165    }
166}
167
168impl Components for PengRobinson {
169    fn components(&self) -> usize {
170        self.parameters.b.len()
171    }
172
173    fn subset(&self, component_list: &[usize]) -> Self {
174        Self::new(Arc::new(self.parameters.subset(component_list)))
175    }
176}
177
178impl Residual for PengRobinson {
179    fn compute_max_density(&self, moles: &Array1<f64>) -> f64 {
180        let b = (moles * &self.parameters.b).sum() / moles.sum();
181        0.9 / b
182    }
183
184    fn residual_helmholtz_energy<D: DualNum<f64> + Copy>(&self, state: &StateHD<D>) -> D {
185        let p = &self.parameters;
186        let x = &state.molefracs;
187        let ak = (&p.tc.mapv(|tc| (D::one() - (state.temperature / tc).sqrt())) * &p.kappa + 1.0)
188            .mapv(|x| x.powi(2))
189            * &p.a;
190
191        // Mixing rules
192        let mut ak_mix = D::zero();
193        for i in 0..ak.len() {
194            for j in 0..ak.len() {
195                ak_mix += (ak[i] * ak[j]).sqrt() * (x[i] * x[j] * (1.0 - p.k_ij[(i, j)]));
196            }
197        }
198        let b = (x * &p.b).sum();
199
200        // Helmholtz energy
201        let n = state.moles.sum();
202        let v = state.volume;
203        n * ((v / (v - b * n)).ln()
204            - ak_mix / (b * SQRT_2 * 2.0 * state.temperature)
205                * ((v + b * n * (1.0 + SQRT_2)) / (v + b * n * (1.0 - SQRT_2))).ln())
206    }
207
208    fn residual_helmholtz_energy_contributions<D: DualNum<f64> + Copy + ScalarOperand>(
209        &self,
210        state: &StateHD<D>,
211    ) -> Vec<(String, D)> {
212        vec![(
213            "Peng Robinson".to_string(),
214            self.residual_helmholtz_energy(state),
215        )]
216    }
217}
218
219impl Molarweight for PengRobinson {
220    fn molar_weight(&self) -> MolarWeight<Array1<f64>> {
221        &self.parameters.molarweight * (GRAM / MOL)
222    }
223}
224
225#[cfg(test)]
226mod tests {
227    use super::*;
228    use crate::state::{Contributions, State};
229    use crate::{EosResult, SolverOptions, Verbosity};
230    use approx::*;
231    use quantity::{KELVIN, PASCAL};
232    use std::sync::Arc;
233
234    fn pure_record_vec() -> Vec<PureRecord<PengRobinsonRecord>> {
235        let records = r#"[
236            {
237                "identifier": {
238                    "cas": "74-98-6",
239                    "name": "propane",
240                    "iupac_name": "propane",
241                    "smiles": "CCC",
242                    "inchi": "InChI=1/C3H8/c1-3-2/h3H2,1-2H3",
243                    "formula": "C3H8"
244                },
245                "model_record": {
246                    "tc": 369.96,
247                    "pc": 4250000.0,
248                    "acentric_factor": 0.153
249                },
250                "molarweight": 44.0962
251            },
252            {
253                "identifier": {
254                    "cas": "106-97-8",
255                    "name": "butane",
256                    "iupac_name": "butane",
257                    "smiles": "CCCC",
258                    "inchi": "InChI=1/C4H10/c1-3-4-2/h3-4H2,1-2H3",
259                    "formula": "C4H10"
260                },
261                "model_record": {
262                    "tc": 425.2,
263                    "pc": 3800000.0,
264                    "acentric_factor": 0.199
265                },
266                "molarweight": 58.123
267            }
268        ]"#;
269        serde_json::from_str(records).expect("Unable to parse json.")
270    }
271
272    #[test]
273    fn peng_robinson() -> EosResult<()> {
274        let mixture = pure_record_vec();
275        let propane = mixture[0].clone();
276        let tc = propane.model_record.tc;
277        let pc = propane.model_record.pc;
278        let parameters = PengRobinsonParameters::new_pure(propane)?;
279        let pr = Arc::new(PengRobinson::new(Arc::new(parameters)));
280        let options = SolverOptions::new().verbosity(Verbosity::Iter);
281        let cp = State::critical_point(&pr, None, None, options)?;
282        println!("{} {}", cp.temperature, cp.pressure(Contributions::Total));
283        assert_relative_eq!(cp.temperature, tc * KELVIN, max_relative = 1e-4);
284        assert_relative_eq!(
285            cp.pressure(Contributions::Total),
286            pc * PASCAL,
287            max_relative = 1e-4
288        );
289        Ok(())
290    }
291}
feos_core/cubic.rs

feos_core/
cubic.rs
