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

feos_core/
cubic.rs
