From ed4d907c542ea159d21f8a9186391449849a393f Mon Sep 17 00:00:00 2001 From: Gernot Bauer Date: Thu, 23 May 2024 17:46:52 +0200 Subject: [PATCH 1/4] started impl for more generic cubic EoSs --- Cargo.toml | 4 +- src/cubic/alpha.rs | 66 ++++++++++++ src/cubic/mixing_rules.rs | 7 ++ src/cubic/mod.rs | 217 ++++++++++++++++++++++++++++++++++++++ src/lib.rs | 2 + 5 files changed, 295 insertions(+), 1 deletion(-) create mode 100644 src/cubic/alpha.rs create mode 100644 src/cubic/mixing_rules.rs create mode 100644 src/cubic/mod.rs diff --git a/Cargo.toml b/Cargo.toml index 2c15ff7ef..3d6d6161c 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -40,6 +40,7 @@ indexmap = "2.0" rayon = { version = "1.7", optional = true } itertools = "0.13" typenum = "1.16" +enum_dispatch = "0.3.13" [dependencies.pyo3] version = "0.21" @@ -73,9 +74,10 @@ uvtheory = ["lazy_static"] pets = [] saftvrqmie = [] saftvrmie = [] +cubic = [] rayon = ["dep:rayon", "ndarray/rayon", "feos-core/rayon", "feos-dft?/rayon"] python = ["pyo3", "numpy", "quantity/python", "feos-core/python", "feos-dft?/python", "rayon"] -all_models = ["dft", "estimator", "pcsaft", "epcsaft", "gc_pcsaft", "uvtheory", "pets", "saftvrqmie", "saftvrmie"] +all_models = ["dft", "estimator", "pcsaft", "epcsaft", "gc_pcsaft", "uvtheory", "pets", "saftvrqmie", "saftvrmie", "cubic"] [[bench]] name = "state_properties" diff --git a/src/cubic/alpha.rs b/src/cubic/alpha.rs new file mode 100644 index 000000000..eae810db3 --- /dev/null +++ b/src/cubic/alpha.rs @@ -0,0 +1,66 @@ +use enum_dispatch::enum_dispatch; +use ndarray::Array1; +use num_dual::DualNum; + +#[enum_dispatch] +trait AlphaFunction { + fn alpha>( + &self, + acentric_factor: &Array1, + reduced_temperature: &Array1, + ) -> Array1; +} + +struct Soave; + +impl AlphaFunction for Soave { + fn alpha>( + &self, + acentric_factor: &Array1, + reduced_temperature: &Array1, + ) -> Array1 { + let m = 0.48 + acentric_factor * (1.574 - acentric_factor * 0.176); + ((-reduced_temperature.mapv(|t| t.sqrt()) + 1.0) * m + 1.0).mapv(|a| a.powi(2)) + } +} + +/// Improved parameterization of the Soave alpha function for Redlich-Kwong equation of state. +/// +/// https://doi.org/10.1016/j.fluid.2018.12.007 +struct SoaveRedlichKwong2019; + +impl AlphaFunction for SoaveRedlichKwong2019 { + fn alpha>( + &self, + acentric_factor: &Array1, + reduced_temperature: &Array1, + ) -> Array1 { + let m = 0.481 + + acentric_factor * (1.5963 - acentric_factor * (0.2963 - acentric_factor * 0.1223)); + ((-reduced_temperature.mapv(|t| t.sqrt()) + 1.0) * m + 1.0).mapv(|a| a.powi(2)) + } +} + +/// Improved parameterization of the Soave alpha function for Peng-Robinson equation of state. +/// +/// https://doi.org/10.1016/j.fluid.2018.12.007 +struct SoavePengRobinson2019; + +impl AlphaFunction for SoavePengRobinson2019 { + fn alpha>( + &self, + acentric_factor: &Array1, + reduced_temperature: &Array1, + ) -> Array1 { + let m = 0.481 + + acentric_factor * (1.5963 - acentric_factor * (0.2963 - acentric_factor * 0.1223)); + ((-reduced_temperature.mapv(|t| t.sqrt()) + 1.0) * m + 1.0).mapv(|a| a.powi(2)) + } +} + +#[enum_dispatch(AlphaFunction)] +enum Alpha { + Soave, + SoaveRedlichKwong2019, + SoavePengRobinson2019, +} diff --git a/src/cubic/mixing_rules.rs b/src/cubic/mixing_rules.rs new file mode 100644 index 000000000..2d74fbfe6 --- /dev/null +++ b/src/cubic/mixing_rules.rs @@ -0,0 +1,7 @@ +use feos_core::parameter::Parameter; +use ndarray::Array1; +use num_dual::DualNum; + +trait MixingRuleFunction { + fn mixing_rule, P: Parameter>(parameters: &P, molefracs: &Array1) -> (D, D); +} diff --git a/src/cubic/mod.rs b/src/cubic/mod.rs new file mode 100644 index 000000000..c93c3a60c --- /dev/null +++ b/src/cubic/mod.rs @@ -0,0 +1,217 @@ +use feos_core::parameter::{Identifier, Parameter, ParameterError, PureRecord}; +use feos_core::si::{MolarWeight, GRAM, MOL}; +use feos_core::{Components, Residual, StateHD}; +use ndarray::{Array1, Array2, ScalarOperand}; +use num_dual::DualNum; +use serde::{Deserialize, Serialize}; +use std::f64::consts::SQRT_2; +use std::fmt; +use std::sync::Arc; + +mod alpha; +mod mixing_rules; + +const KB_A3: f64 = 13806490.0; + +/// Cubic parameters for a single substance. +#[derive(Serialize, Deserialize, Debug, Clone, Default)] +pub struct CubicRecord { + /// critical temperature in Kelvin + tc: f64, + /// critical pressure in Pascal + pc: f64, + /// acentric factor + acentric_factor: f64, +} + +impl CubicRecord { + /// Create a new pure substance record for the Cubic equation of state. + pub fn new(tc: f64, pc: f64, acentric_factor: f64) -> Self { + Self { + tc, + pc, + acentric_factor, + } + } +} + +impl std::fmt::Display for CubicRecord { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + write!(f, "CubicRecord(tc={} K", self.tc)?; + write!(f, ", pc={} Pa", self.pc)?; + write!(f, ", acentric factor={}", self.acentric_factor) + } +} + +/// Cubic parameters for one ore more substances. +pub struct CubicParameters { + /// Critical temperature in Kelvin + tc: Array1, + a: Array1, + b: Array1, + /// Binary interaction parameter + k_ij: Array2, + kappa: Array1, + /// Molar weight in units of g/mol + molarweight: Array1, + /// List of pure component records + pure_records: Vec>, +} + +impl std::fmt::Display for CubicParameters { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + self.pure_records + .iter() + .try_for_each(|pr| writeln!(f, "{}", pr))?; + writeln!(f, "\nk_ij:\n{}", self.k_ij) + } +} + +impl CubicParameters { + /// Build a simple parameter set without binary interaction parameters. + pub fn new_simple( + tc: &[f64], + pc: &[f64], + acentric_factor: &[f64], + molarweight: &[f64], + ) -> Result { + if [pc.len(), acentric_factor.len(), molarweight.len()] + .iter() + .any(|&l| l != tc.len()) + { + return Err(ParameterError::IncompatibleParameters(String::from( + "each component has to have parameters.", + ))); + } + let records = (0..tc.len()) + .map(|i| { + let record = CubicRecord { + tc: tc[i], + pc: pc[i], + acentric_factor: acentric_factor[i], + }; + let id = Identifier::default(); + PureRecord::new(id, molarweight[i], record) + }) + .collect(); + CubicParameters::from_records(records, None) + } +} + +impl Parameter for CubicParameters { + type Pure = CubicRecord; + type Binary = f64; + + /// Creates parameters from pure component records. + fn from_records( + pure_records: Vec>, + binary_records: Option>, + ) -> Result { + let n = pure_records.len(); + + let mut tc = Array1::zeros(n); + let mut a = Array1::zeros(n); + let mut b = Array1::zeros(n); + let mut molarweight = Array1::zeros(n); + let mut kappa = Array1::zeros(n); + + for (i, record) in pure_records.iter().enumerate() { + molarweight[i] = record.molarweight; + let r = &record.model_record; + tc[i] = r.tc; + a[i] = 0.45724 * r.tc.powi(2) * KB_A3 / r.pc; + b[i] = 0.07780 * r.tc * KB_A3 / r.pc; + kappa[i] = 0.37464 + (1.54226 - 0.26992 * r.acentric_factor) * r.acentric_factor; + } + + let k_ij = binary_records.unwrap_or_else(|| Array2::zeros([n; 2])); + + Ok(Self { + tc, + a, + b, + k_ij, + kappa, + molarweight, + pure_records, + }) + } + + fn records(&self) -> (&[PureRecord], Option<&Array2>) { + (&self.pure_records, Some(&self.k_ij)) + } +} + +/// A simple version of the Cubic equation of state. +pub struct Cubic { + /// Parameters + parameters: Arc, +} + +impl Cubic { + /// Create a new equation of state from a set of parameters. + pub fn new(parameters: Arc) -> Self { + Self { parameters } + } +} + +impl fmt::Display for Cubic { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write!(f, "Peng Robinson") + } +} + +impl Components for Cubic { + fn components(&self) -> usize { + self.parameters.b.len() + } + + fn subset(&self, component_list: &[usize]) -> Self { + Self::new(Arc::new(self.parameters.subset(component_list))) + } +} + +impl Residual for Cubic { + fn compute_max_density(&self, moles: &Array1) -> f64 { + let b = (moles * &self.parameters.b).sum() / moles.sum(); + 0.9 / b + } + + fn residual_helmholtz_energy + Copy>(&self, state: &StateHD) -> D { + let p = &self.parameters; + let x = &state.molefracs; + let ak = (&p.tc.mapv(|tc| (D::one() - (state.temperature / tc).sqrt())) * &p.kappa + 1.0) + .mapv(|x| x.powi(2)) + * &p.a; + + // Mixing rules + let mut ak_mix = D::zero(); + for i in 0..ak.len() { + for j in 0..ak.len() { + ak_mix += (ak[i] * ak[j]).sqrt() * (x[i] * x[j] * (1.0 - p.k_ij[(i, j)])); + } + } + let b = (x * &p.b).sum(); + + // Helmholtz energy + let n = state.moles.sum(); + let v = state.volume; + n * ((v / (v - b * n)).ln() + - ak_mix / (b * SQRT_2 * 2.0 * state.temperature) + * ((v + b * n * (1.0 + SQRT_2)) / (v + b * n * (1.0 - SQRT_2))).ln()) + } + + fn residual_helmholtz_energy_contributions + Copy + ScalarOperand>( + &self, + state: &StateHD, + ) -> Vec<(String, D)> { + vec![( + "Peng Robinson".to_string(), + self.residual_helmholtz_energy(state), + )] + } + + fn molar_weight(&self) -> MolarWeight> { + &self.parameters.molarweight * (GRAM / MOL) + } +} diff --git a/src/lib.rs b/src/lib.rs index 2f5659b2d..4f473a122 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -50,6 +50,8 @@ pub mod association; pub mod hard_sphere; // models +#[cfg(feature = "cubic")] +pub mod cubic; #[cfg(feature = "epcsaft")] pub mod epcsaft; #[cfg(feature = "gc_pcsaft")] From cd5cfbe5508e5d19965e657e62badfdb9d408207 Mon Sep 17 00:00:00 2001 From: Gernot Bauer Date: Fri, 7 Jun 2024 15:08:46 +0200 Subject: [PATCH 2/4] added several alpha functions --- feos-core/src/cubic.rs | 2 + src/cubic/alpha.rs | 66 ------ src/cubic/alpha/mathias_copeman.rs | 52 +++++ src/cubic/alpha/mod.rs | 47 ++++ src/cubic/alpha/soave.rs | 173 ++++++++++++++ src/cubic/alpha/twu.rs | 127 +++++++++++ src/cubic/mixing_rules.rs | 64 +++++- src/cubic/mod.rs | 349 +++++++++++++++++------------ src/cubic/parameters.rs | 217 ++++++++++++++++++ src/cubic/python.rs | 293 ++++++++++++++++++++++++ src/eos.rs | 4 + src/python/cubic.rs | 2 +- src/python/eos.rs | 91 +++++++- src/python/mod.rs | 7 +- 14 files changed, 1275 insertions(+), 219 deletions(-) delete mode 100644 src/cubic/alpha.rs create mode 100644 src/cubic/alpha/mathias_copeman.rs create mode 100644 src/cubic/alpha/mod.rs create mode 100644 src/cubic/alpha/soave.rs create mode 100644 src/cubic/alpha/twu.rs create mode 100644 src/cubic/parameters.rs create mode 100644 src/cubic/python.rs diff --git a/feos-core/src/cubic.rs b/feos-core/src/cubic.rs index 1b178ef51..bb51100a3 100644 --- a/feos-core/src/cubic.rs +++ b/feos-core/src/cubic.rs @@ -200,6 +200,8 @@ impl Residual for PengRobinson { // Helmholtz energy let n = state.moles.sum(); let v = state.volume; + dbg!(&b); + dbg!(&ak_mix); n * ((v / (v - b * n)).ln() - ak_mix / (b * SQRT_2 * 2.0 * state.temperature) * ((v + b * n * (1.0 + SQRT_2)) / (v + b * n * (1.0 - SQRT_2))).ln()) diff --git a/src/cubic/alpha.rs b/src/cubic/alpha.rs deleted file mode 100644 index eae810db3..000000000 --- a/src/cubic/alpha.rs +++ /dev/null @@ -1,66 +0,0 @@ -use enum_dispatch::enum_dispatch; -use ndarray::Array1; -use num_dual::DualNum; - -#[enum_dispatch] -trait AlphaFunction { - fn alpha>( - &self, - acentric_factor: &Array1, - reduced_temperature: &Array1, - ) -> Array1; -} - -struct Soave; - -impl AlphaFunction for Soave { - fn alpha>( - &self, - acentric_factor: &Array1, - reduced_temperature: &Array1, - ) -> Array1 { - let m = 0.48 + acentric_factor * (1.574 - acentric_factor * 0.176); - ((-reduced_temperature.mapv(|t| t.sqrt()) + 1.0) * m + 1.0).mapv(|a| a.powi(2)) - } -} - -/// Improved parameterization of the Soave alpha function for Redlich-Kwong equation of state. -/// -/// https://doi.org/10.1016/j.fluid.2018.12.007 -struct SoaveRedlichKwong2019; - -impl AlphaFunction for SoaveRedlichKwong2019 { - fn alpha>( - &self, - acentric_factor: &Array1, - reduced_temperature: &Array1, - ) -> Array1 { - let m = 0.481 - + acentric_factor * (1.5963 - acentric_factor * (0.2963 - acentric_factor * 0.1223)); - ((-reduced_temperature.mapv(|t| t.sqrt()) + 1.0) * m + 1.0).mapv(|a| a.powi(2)) - } -} - -/// Improved parameterization of the Soave alpha function for Peng-Robinson equation of state. -/// -/// https://doi.org/10.1016/j.fluid.2018.12.007 -struct SoavePengRobinson2019; - -impl AlphaFunction for SoavePengRobinson2019 { - fn alpha>( - &self, - acentric_factor: &Array1, - reduced_temperature: &Array1, - ) -> Array1 { - let m = 0.481 - + acentric_factor * (1.5963 - acentric_factor * (0.2963 - acentric_factor * 0.1223)); - ((-reduced_temperature.mapv(|t| t.sqrt()) + 1.0) * m + 1.0).mapv(|a| a.powi(2)) - } -} - -#[enum_dispatch(AlphaFunction)] -enum Alpha { - Soave, - SoaveRedlichKwong2019, - SoavePengRobinson2019, -} diff --git a/src/cubic/alpha/mathias_copeman.rs b/src/cubic/alpha/mathias_copeman.rs new file mode 100644 index 000000000..2b1f3cd3b --- /dev/null +++ b/src/cubic/alpha/mathias_copeman.rs @@ -0,0 +1,52 @@ +use std::sync::Arc; + +use feos_core::{parameter::ParameterError, EosResult}; +use ndarray::{Array1, Zip}; +use num_dual::DualNum; +use serde::{Deserialize, Serialize}; + +use crate::cubic::parameters::CubicParameters; + +use super::AlphaFunction; + +#[derive(Serialize, Deserialize, Clone, Debug)] +pub struct MathiasCopeman(pub Vec<[f64; 3]>); + +impl AlphaFunction for MathiasCopeman { + #[inline] + fn alpha + Copy>( + &self, + _: &Array1, + reduced_temperature: &Array1, + ) -> Array1 { + Zip::from(reduced_temperature) + .and(&self.0) + .map_collect(|tr, c| { + let trsq = -tr.sqrt() + 1.0; + let a1 = trsq + c[0] + 1.0; + let a2 = match tr { + tr if tr.re() < 1.0 => trsq * (c[1] + c[2]), + _ => D::zero(), + }; + (a1 + a2).powi(2) + }) + } + + fn validate(&self, parameters: &Arc) -> EosResult<()> { + if self.0.len() == parameters.tc.len() { + Ok(()) + } else { + Err(ParameterError::IncompatibleParameters( + format!( + "Mathias Copeman alpha function was initialized for {} components, but the equation of state contains {}.", + self.0.len(), parameters.tc.len() + ) + ).into()) + } + } + + fn subset(&self, component_list: &[usize]) -> Self { + let mi = component_list.iter().map(|&i| self.0[i]).collect(); + Self(mi) + } +} diff --git a/src/cubic/alpha/mod.rs b/src/cubic/alpha/mod.rs new file mode 100644 index 000000000..3ec1bafc1 --- /dev/null +++ b/src/cubic/alpha/mod.rs @@ -0,0 +1,47 @@ +use std::sync::Arc; + +use enum_dispatch::enum_dispatch; +use feos_core::EosResult; +use ndarray::{Array1, ScalarOperand}; +use num_dual::DualNum; + +use super::parameters::CubicParameters; + +mod soave; +pub use soave::{ + PengRobinson1976, PengRobinson1978, PengRobinson2019, RedlichKwong1972, RedlichKwong2019, Soave, +}; +mod mathias_copeman; +pub use mathias_copeman::MathiasCopeman; +mod twu; +pub use twu::{GeneralizedTwu, Twu}; + +#[enum_dispatch] +pub trait AlphaFunction { + fn alpha + Copy + ScalarOperand>( + &self, + acentric_factor: &Array1, + reduced_temperature: &Array1, + ) -> Array1; + + /// Check for validity of alpha function against parameters, e.g. + /// to assert that the number of components match. + fn validate(&self, parameters: &Arc) -> EosResult<()>; + + /// Generate the alpha function for a subset of components. + fn subset(&self, component_list: &[usize]) -> Self; +} + +#[enum_dispatch(AlphaFunction)] +#[derive(Debug, Clone)] +pub enum Alpha { + Soave, + PengRobinson1976, + PengRobinson1978, + PengRobinson2019, + RedlichKwong1972, + RedlichKwong2019, + MathiasCopeman, + GeneralizedTwu, + Twu, +} diff --git a/src/cubic/alpha/soave.rs b/src/cubic/alpha/soave.rs new file mode 100644 index 000000000..b8de68454 --- /dev/null +++ b/src/cubic/alpha/soave.rs @@ -0,0 +1,173 @@ +use std::sync::Arc; + +use feos_core::EosResult; +use ndarray::{Array1, Zip}; +use num_dual::DualNum; +use serde::{Deserialize, Serialize}; + +use crate::cubic::parameters::CubicParameters; + +use super::AlphaFunction; + +#[derive(Serialize, Deserialize, Clone, Debug)] +/// Generic version of Soave's function using 3rd order polynomial. +pub struct Soave { + /// coefficients for m-polynomial + mi: Array1, +} + +impl Soave { + pub fn new(mi: Array1) -> Self { + Soave { mi } + } +} + +impl AlphaFunction for Soave { + #[inline] + fn alpha>( + &self, + acentric_factor: &Array1, + reduced_temperature: &Array1, + ) -> Array1 { + let m = self + .mi + .iter() + .enumerate() + .fold(Array1::zeros(acentric_factor.len()), |m, (i, &mi)| { + &m + acentric_factor.mapv(|w| w.powi(i as i32)) * mi + }); + ((-reduced_temperature.mapv(|t| t.sqrt()) + 1.0) * m + 1.0).mapv(|a| a.powi(2)) + } + + fn validate(&self, _: &Arc) -> EosResult<()> { + Ok(()) + } + + fn subset(&self, _: &[usize]) -> Self { + self.clone() + } +} + +#[derive(Serialize, Deserialize, Clone, Debug)] +pub struct RedlichKwong1972; + +impl AlphaFunction for RedlichKwong1972 { + #[inline] + fn alpha>( + &self, + acentric_factor: &Array1, + reduced_temperature: &Array1, + ) -> Array1 { + let m = acentric_factor.mapv(|w| 0.48 + w * (1.574 - w * 0.176)); + ((-reduced_temperature.mapv(|t| t.sqrt()) + 1.0) * m + 1.0).mapv(|a| a.powi(2)) + } + + fn validate(&self, _: &Arc) -> EosResult<()> { + Ok(()) + } + + fn subset(&self, _: &[usize]) -> Self { + Self + } +} + +#[derive(Serialize, Deserialize, Clone, Debug)] +pub struct PengRobinson1976; + +impl AlphaFunction for PengRobinson1976 { + #[inline] + fn alpha>( + &self, + acentric_factor: &Array1, + reduced_temperature: &Array1, + ) -> Array1 { + let m = acentric_factor.mapv(|w| 0.37464 + w * (1.54226 - w * 0.26992)); + ((-reduced_temperature.mapv(|t| t.sqrt()) + 1.0) * m + 1.0).mapv(|a| a.powi(2)) + } + + fn validate(&self, _: &Arc) -> EosResult<()> { + Ok(()) + } + + fn subset(&self, _: &[usize]) -> Self { + Self + } +} + +#[derive(Serialize, Deserialize, Clone, Debug)] +pub struct PengRobinson1978; + +impl AlphaFunction for PengRobinson1978 { + #[inline] + fn alpha + Copy>( + &self, + acentric_factor: &Array1, + reduced_temperature: &Array1, + ) -> Array1 { + Zip::from(acentric_factor) + .and(reduced_temperature) + .map_collect(|&w, &tr| { + let m = if w <= 0.491 { + 0.37464 + w * (1.54226 - w * 0.26992) + } else { + // use higher-order polynomial if w > w(n-decane) + 0.379642 + w * (1.48503 + w * (-0.164423 + w * 0.016666)) + }; + ((-tr.sqrt() + 1.0) * m + 1.0).powi(2) + }) + } + + fn validate(&self, _: &Arc) -> EosResult<()> { + Ok(()) + } + + fn subset(&self, _: &[usize]) -> Self { + Self + } +} + +#[derive(Serialize, Deserialize, Clone, Debug)] +pub struct RedlichKwong2019; + +impl AlphaFunction for RedlichKwong2019 { + #[inline] + fn alpha>( + &self, + acentric_factor: &Array1, + reduced_temperature: &Array1, + ) -> Array1 { + let m = acentric_factor.mapv(|w| 0.481 + w * (1.5963 + w * (-0.2963 + w * 0.1223))); + ((-reduced_temperature.mapv(|t| t.sqrt()) + 1.0) * m + 1.0).mapv(|a| a.powi(2)) + } + + fn validate(&self, _: &Arc) -> EosResult<()> { + Ok(()) + } + + fn subset(&self, _: &[usize]) -> Self { + Self + } +} + +#[derive(Serialize, Deserialize, Clone, Debug)] +pub struct PengRobinson2019; + +impl AlphaFunction for PengRobinson2019 { + #[inline] + fn alpha>( + &self, + acentric_factor: &Array1, + reduced_temperature: &Array1, + ) -> Array1 { + let m = acentric_factor.mapv(|w| 0.3919 + w * (1.4996 + w * (-0.2721 + w * 0.1063))); + ((-reduced_temperature.mapv(|t| t.sqrt()) + 1.0) * m + 1.0).mapv(|a| a.powi(2)) + } + + fn validate(&self, _: &Arc) -> EosResult<()> { + Ok(()) + } + + fn subset(&self, _: &[usize]) -> Self { + Self + } +} diff --git a/src/cubic/alpha/twu.rs b/src/cubic/alpha/twu.rs new file mode 100644 index 000000000..d30083292 --- /dev/null +++ b/src/cubic/alpha/twu.rs @@ -0,0 +1,127 @@ +use std::sync::Arc; + +use feos_core::{parameter::ParameterError, EosResult}; +use itertools::izip; +use ndarray::{Array1, Zip}; +use num_dual::DualNum; +use serde::{Deserialize, Serialize}; + +use crate::cubic::parameters::CubicParameters; + +use super::AlphaFunction; + +/// Generalized version of the Twu alpha function (1995). +/// +/// Different parameters are used for sub- and supercritical conditions. +#[derive(Serialize, Deserialize, Clone, Debug)] +pub struct GeneralizedTwu([[[f64; 3]; 2]; 2]); + +impl GeneralizedTwu { + pub fn redlich_kwong() -> Self { + GeneralizedTwu([ + [ + [2.496441 * (0.919422 - 1.0), 0.141599, 2.496441 * 0.919422], + [3.291790 * (0.799457 - 1.0), 0.500315, 3.291790 * 0.799457], + ], + [ + [-0.2 * (6.500018 - 1.0), 0.441411, -0.2 * 6.500018], + [-8.0 * (1.289098 - 1.0), 0.032580, -8.0 * 1.289098], + ], + ]) + } + + pub fn peng_robinson() -> Self { + GeneralizedTwu([ + [ + [1.948150 * (0.911807 - 1.0), 0.125283, 1.948150 * 0.911807], + [2.812520 * (0.784054 - 1.0), 0.511614, 2.812520 * 0.784054], + ], + [ + [-0.2 * (4.963070 - 1.0), 0.401219, -0.2 * 4.963070], + [-0.8 * (1.248089 - 1.0), 0.024955, -0.8 * 1.248089], + ], + ]) + } +} + +impl AlphaFunction for GeneralizedTwu { + #[inline] + fn alpha + Copy>( + &self, + acentric_factor: &Array1, + reduced_temperature: &Array1, + ) -> Array1 { + Zip::from(acentric_factor) + .and(reduced_temperature) + .map_collect(|&w, &tr| { + if tr.re() <= 1.0 { + let [nm_m1_0, l0, nm0] = self.0[0][0]; + let [nm_m1_1, l1, nm1] = self.0[0][1]; + let a0 = tr.powf(nm_m1_0) * ((-tr.powf(nm0) + 1.0) * l0).exp(); + let a1 = tr.powf(nm_m1_1) * ((-tr.powf(nm1) + 1.0) * l1).exp(); + a0 + (a1 - a0) * w + } else { + let [nm_m1_0, l0, nm0] = self.0[1][0]; + let [nm_m1_1, l1, nm1] = self.0[1][1]; + let a0 = tr.powf(nm_m1_0) * ((-tr.powf(nm0) + 1.0) * l0).exp(); + let a1 = tr.powf(nm_m1_1) * ((-tr.powf(nm1) + 1.0) * l1).exp(); + a0 + (a1 - a0) * w + } + }) + } + + fn validate(&self, _: &Arc) -> EosResult<()> { + Ok(()) + } + + fn subset(&self, _: &[usize]) -> Self { + self.clone() + } +} + +/// Generalized version of the Twu alpha function (1995). +/// +/// Different parameters are used for sub- and supercritical conditions. +#[derive(Serialize, Deserialize, Clone, Debug)] +pub struct Twu(Vec<[f64; 3]>); + +impl Twu { + pub fn new(l: Vec, m: Vec, n: Option>) -> Self { + let _n = n.unwrap_or(vec![2.0; l.len()]); + let input = izip!(l, m, _n) + .map(|(l, m, n)| [n * (m - 1.0), l, n * m]) + .collect(); + Self(input) + } +} + +impl AlphaFunction for Twu { + #[inline] + fn alpha + Copy>( + &self, + _: &Array1, + reduced_temperature: &Array1, + ) -> Array1 { + izip!(reduced_temperature, &self.0) + .map(|(tr, &[nmm1, l, nm])| tr.powf(nmm1) * ((-tr.powf(nm) + 1.0) * l).exp()) + .collect() + } + + fn validate(&self, parameters: &Arc) -> EosResult<()> { + if self.0.len() == parameters.tc.len() { + Ok(()) + } else { + Err(ParameterError::IncompatibleParameters( + format!( + "Twu alpha function was initialized for {} components, but the equation of state contains {}.", + self.0.len(), parameters.tc.len() + ) + ).into()) + } + } + + fn subset(&self, component_list: &[usize]) -> Self { + let c = component_list.iter().map(|&i| self.0[i]).collect(); + Self(c) + } +} diff --git a/src/cubic/mixing_rules.rs b/src/cubic/mixing_rules.rs index 2d74fbfe6..8bfdc356d 100644 --- a/src/cubic/mixing_rules.rs +++ b/src/cubic/mixing_rules.rs @@ -1,7 +1,63 @@ -use feos_core::parameter::Parameter; -use ndarray::Array1; +use enum_dispatch::enum_dispatch; +use feos_core::StateHD; +use ndarray::ScalarOperand; use num_dual::DualNum; -trait MixingRuleFunction { - fn mixing_rule, P: Parameter>(parameters: &P, molefracs: &Array1) -> (D, D); +use super::{alpha::AlphaFunction, Cubic}; + +/// Parameters of cubics +pub struct MixtureParameters { + /// attractive contribution divided by R + pub a: D, + /// repulsive contribution + pub b: D, + /// volume translation + pub c: D, +} + +#[enum_dispatch] +pub trait MixingRuleFunction { + fn apply + Copy + ScalarOperand>( + &self, + cubic: &Cubic, + state: &StateHD, + ) -> MixtureParameters; +} + +/// Quadratic summation over a and b. +#[derive(Debug, Clone)] +pub struct Quadratic; + +impl MixingRuleFunction for Quadratic { + fn apply + Copy + ScalarOperand>( + &self, + cubic: &Cubic, + state: &StateHD, + ) -> MixtureParameters { + let p = &cubic.parameters; + let pc = &cubic.critical_parameters; + let n = p.tc.len(); + let tr = p.tc.mapv(|tc| state.temperature / tc); + let at = cubic.options.alpha.alpha(&p.acentric_factor, &tr) * &pc.ac; + let mut a = D::zero(); + let mut b = D::zero(); + for i in 0..n { + let xi = state.molefracs[i]; + let ai = at[i]; + let bi = pc.bc[i]; + a += xi * xi * ai; + b += xi * xi * bi; + for j in i + 1..n { + a += xi * state.molefracs[j] * (ai * at[j]).sqrt() * (1.0 - p.k_ij[[i, j]]) * 2.0; + b += xi * state.molefracs[j] * (bi + pc.bc[j]) * 0.5 * (1.0 - p.l_ij[[i, j]]) * 2.0; + } + } + MixtureParameters { a, b, c: D::zero() } + } +} + +#[enum_dispatch(MixingRuleFunction)] +#[derive(Debug, Clone)] +pub enum MixingRule { + Quadratic, } diff --git a/src/cubic/mod.rs b/src/cubic/mod.rs index c93c3a60c..6f33bdc63 100644 --- a/src/cubic/mod.rs +++ b/src/cubic/mod.rs @@ -1,217 +1,274 @@ -use feos_core::parameter::{Identifier, Parameter, ParameterError, PureRecord}; +use feos_core::parameter::Parameter; use feos_core::si::{MolarWeight, GRAM, MOL}; -use feos_core::{Components, Residual, StateHD}; -use ndarray::{Array1, Array2, ScalarOperand}; +use feos_core::{Components, EosResult, Residual, StateHD}; +use mixing_rules::{MixingRule, MixingRuleFunction, MixtureParameters, Quadratic}; +use ndarray::{Array1, ScalarOperand, Zip}; use num_dual::DualNum; -use serde::{Deserialize, Serialize}; use std::f64::consts::SQRT_2; use std::fmt; use std::sync::Arc; mod alpha; +use alpha::{Alpha, AlphaFunction, PengRobinson1976, RedlichKwong1972}; mod mixing_rules; +mod parameters; +use parameters::CubicParameters; +// mod volume_translation; + +#[cfg(feature = "python")] +pub(crate) mod python; const KB_A3: f64 = 13806490.0; -/// Cubic parameters for a single substance. -#[derive(Serialize, Deserialize, Debug, Clone, Default)] -pub struct CubicRecord { - /// critical temperature in Kelvin - tc: f64, - /// critical pressure in Pascal - pc: f64, - /// acentric factor - acentric_factor: f64, +#[derive(Debug, Clone, Copy)] +pub struct Delta { + d1: f64, + d2: f64, + d12: f64, } -impl CubicRecord { - /// Create a new pure substance record for the Cubic equation of state. - pub fn new(tc: f64, pc: f64, acentric_factor: f64) -> Self { - Self { - tc, - pc, - acentric_factor, +impl From<(f64, f64)> for Delta { + fn from(value: (f64, f64)) -> Self { + Delta { + d1: value.0, + d2: value.1, + d12: value.0 - value.1, } } } -impl std::fmt::Display for CubicRecord { - fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { - write!(f, "CubicRecord(tc={} K", self.tc)?; - write!(f, ", pc={} Pa", self.pc)?; - write!(f, ", acentric factor={}", self.acentric_factor) +impl Delta { + // Calculate universal critical constants from universal cubic parameters. + // + // See https://doi.org/10.1016/j.fluid.2012.05.008 + fn critical_constants(&self) -> (f64, f64) { + let (r1, r2) = (-self.d1, -self.d2); + let eta_c = 1.0 + / (((1.0 - r1) * (1.0 - r2).powi(2)).cbrt() + + ((1.0 - r2) * (1.0 - r1).powi(2)).cbrt() + + 1.0); + let omega_a = (1.0 - eta_c * r1) * (1.0 - eta_c * r2) / (1.0 - eta_c) + * (2.0 - eta_c * (r1 + r2)) + / (3.0 - eta_c * (1.0 + r1 + r2)).powi(2); + let omega_b = eta_c / (3.0 - eta_c * (1.0 + r1 + r2)); + (omega_a, omega_b) } } -/// Cubic parameters for one ore more substances. -pub struct CubicParameters { - /// Critical temperature in Kelvin - tc: Array1, - a: Array1, - b: Array1, - /// Binary interaction parameter - k_ij: Array2, - kappa: Array1, - /// Molar weight in units of g/mol - molarweight: Array1, - /// List of pure component records - pure_records: Vec>, +/// Parameters processed using model constants and substance critial data. +#[derive(Debug)] +pub struct CriticalParameters { + ac: Array1, + bc: Array1, + omega_a: f64, + omega_b: f64, } -impl std::fmt::Display for CubicParameters { - fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { - self.pure_records - .iter() - .try_for_each(|pr| writeln!(f, "{}", pr))?; - writeln!(f, "\nk_ij:\n{}", self.k_ij) +impl CriticalParameters { + fn new(p: &Arc, delta: &Delta) -> Self { + let (omega_a, omega_b) = delta.critical_constants(); + let ac = omega_a * &p.tc.mapv(|tc| tc.powi(2)) * KB_A3 / &p.pc; + let bc = omega_b * &p.tc * KB_A3 / &p.pc; + Self { + ac, + bc, + omega_a, + omega_b, + } } -} -impl CubicParameters { - /// Build a simple parameter set without binary interaction parameters. - pub fn new_simple( - tc: &[f64], - pc: &[f64], - acentric_factor: &[f64], - molarweight: &[f64], - ) -> Result { - if [pc.len(), acentric_factor.len(), molarweight.len()] - .iter() - .any(|&l| l != tc.len()) - { - return Err(ParameterError::IncompatibleParameters(String::from( - "each component has to have parameters.", - ))); + fn subset(&self, component_list: &[usize]) -> Self { + let n = component_list.len(); + let mut ac = Array1::zeros(n); + let mut bc = Array1::zeros(n); + Zip::from(&mut ac) + .and(&mut bc) + .and(component_list) + .for_each(|a, b, &i| { + *a = self.ac[i]; + *b = self.bc[i]; + }); + Self { + ac, + bc, + omega_a: self.omega_a, + omega_b: self.omega_b, } - let records = (0..tc.len()) - .map(|i| { - let record = CubicRecord { - tc: tc[i], - pc: pc[i], - acentric_factor: acentric_factor[i], - }; - let id = Identifier::default(); - PureRecord::new(id, molarweight[i], record) - }) - .collect(); - CubicParameters::from_records(records, None) } } -impl Parameter for CubicParameters { - type Pure = CubicRecord; - type Binary = f64; - - /// Creates parameters from pure component records. - fn from_records( - pure_records: Vec>, - binary_records: Option>, - ) -> Result { - let n = pure_records.len(); - - let mut tc = Array1::zeros(n); - let mut a = Array1::zeros(n); - let mut b = Array1::zeros(n); - let mut molarweight = Array1::zeros(n); - let mut kappa = Array1::zeros(n); - - for (i, record) in pure_records.iter().enumerate() { - molarweight[i] = record.molarweight; - let r = &record.model_record; - tc[i] = r.tc; - a[i] = 0.45724 * r.tc.powi(2) * KB_A3 / r.pc; - b[i] = 0.07780 * r.tc * KB_A3 / r.pc; - kappa[i] = 0.37464 + (1.54226 - 0.26992 * r.acentric_factor) * r.acentric_factor; - } +#[derive(Debug, Clone)] +pub struct CubicOptions { + pub(crate) alpha: Alpha, + pub(crate) mixing: MixingRule, + pub(crate) delta: Delta, +} - let k_ij = binary_records.unwrap_or_else(|| Array2::zeros([n; 2])); +/// A generic cubic equation of state. +pub struct Cubic { + /// Parameters + pub parameters: Arc, + pub options: CubicOptions, + /// processed parameters using model and substance critical data + pub critical_parameters: CriticalParameters, +} +impl Cubic { + /// Generic cubic equation of state with adjustable universal constants. + pub fn new(parameters: Arc, options: CubicOptions) -> EosResult { + let p = CriticalParameters::new(¶meters, &options.delta); + options.alpha.validate(¶meters)?; Ok(Self { - tc, - a, - b, - k_ij, - kappa, - molarweight, - pure_records, + parameters, + options, + critical_parameters: p, }) } - fn records(&self) -> (&[PureRecord], Option<&Array2>) { - (&self.pure_records, Some(&self.k_ij)) + /// Peng Robinson equation of state. + /// + /// Universal constants: + /// - $\delta_1 = 1 + \sqrt{2}$ + /// - $\delta_2 = 1 - \sqrt{2}$ + /// + /// If no options are supplied, the following is used: + /// - alpha function: Peng Robinson (1976) + /// - mixing rules: quadratic mixing + pub fn peng_robinson( + parameters: Arc, + alpha: Option, + mixing: Option, + ) -> EosResult { + let delta: Delta = (1.0 + SQRT_2, 1.0 - SQRT_2).into(); + let p = CriticalParameters::new(¶meters, &delta); + let options = CubicOptions { + alpha: alpha.unwrap_or(PengRobinson1976.into()), + mixing: mixing.unwrap_or(Quadratic.into()), + delta, + }; + options.alpha.validate(¶meters)?; + Ok(Self { + parameters, + options, + critical_parameters: p, + }) } -} - -/// A simple version of the Cubic equation of state. -pub struct Cubic { - /// Parameters - parameters: Arc, -} -impl Cubic { - /// Create a new equation of state from a set of parameters. - pub fn new(parameters: Arc) -> Self { - Self { parameters } + /// Create equation of state of (Suave) Redlich Kwong. + /// + /// Universal constants: + /// - $\delta_1 = 1$ + /// - $\delta_2 = 0$ + /// + /// If no options are supplied, the following is used: + /// - alpha function: Soave (1972) + /// - mixing rules: quadratic mixing + pub fn redlich_kwong( + parameters: Arc, + alpha: Option, + mixing: Option, + ) -> EosResult { + let delta: Delta = (1.0, 0.0).into(); + let p = CriticalParameters::new(¶meters, &delta); + let options = CubicOptions { + alpha: alpha.unwrap_or(RedlichKwong1972.into()), + mixing: mixing.unwrap_or(Quadratic.into()), + delta, + }; + options.alpha.validate(¶meters)?; + Ok(Self { + parameters, + options, + critical_parameters: p, + }) } } impl fmt::Display for Cubic { fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!(f, "Peng Robinson") + write!(f, "cubic") } } impl Components for Cubic { fn components(&self) -> usize { - self.parameters.b.len() + self.parameters.tc.len() } fn subset(&self, component_list: &[usize]) -> Self { - Self::new(Arc::new(self.parameters.subset(component_list))) + Self { + parameters: Arc::new(self.parameters.subset(component_list)), + options: self.options.clone(), + critical_parameters: self.critical_parameters.subset(component_list), + } } } impl Residual for Cubic { fn compute_max_density(&self, moles: &Array1) -> f64 { - let b = (moles * &self.parameters.b).sum() / moles.sum(); + let b = (moles * &self.critical_parameters.bc).sum() / moles.sum(); 0.9 / b } - fn residual_helmholtz_energy + Copy>(&self, state: &StateHD) -> D { - let p = &self.parameters; - let x = &state.molefracs; - let ak = (&p.tc.mapv(|tc| (D::one() - (state.temperature / tc).sqrt())) * &p.kappa + 1.0) - .mapv(|x| x.powi(2)) - * &p.a; - - // Mixing rules - let mut ak_mix = D::zero(); - for i in 0..ak.len() { - for j in 0..ak.len() { - ak_mix += (ak[i] * ak[j]).sqrt() * (x[i] * x[j] * (1.0 - p.k_ij[(i, j)])); - } - } - let b = (x * &p.b).sum(); - - // Helmholtz energy + fn residual_helmholtz_energy + Copy + ScalarOperand>( + &self, + state: &StateHD, + ) -> D { + let MixtureParameters { a, b, c: _ } = self.options.mixing.apply(self, state); let n = state.moles.sum(); let v = state.volume; - n * ((v / (v - b * n)).ln() - - ak_mix / (b * SQRT_2 * 2.0 * state.temperature) - * ((v + b * n * (1.0 + SQRT_2)) / (v + b * n * (1.0 - SQRT_2))).ln()) + let bn = b * n; + n * ((v / (v - bn)).ln() + - a / (b * self.options.delta.d12 * state.temperature) + * ((v + bn * self.options.delta.d1) / (v + bn * self.options.delta.d2)).ln()) } fn residual_helmholtz_energy_contributions + Copy + ScalarOperand>( &self, state: &StateHD, ) -> Vec<(String, D)> { - vec![( - "Peng Robinson".to_string(), - self.residual_helmholtz_energy(state), - )] + vec![("cubic".to_string(), self.residual_helmholtz_energy(state))] } fn molar_weight(&self) -> MolarWeight> { &self.parameters.molarweight * (GRAM / MOL) } } + +#[cfg(test)] +mod tests { + use feos_core::{ + cubic::{PengRobinson, PengRobinsonParameters, PengRobinsonRecord}, + parameter::{Identifier, PureRecord}, + }; + use ndarray::arr1; + use parameters::CubicRecord; + + use super::*; + + #[test] + fn a_res() { + let propane = PureRecord::new( + Identifier::new(None, Some("propane"), None, None, None, None), + 44.0962, + CubicRecord::new(369.96, 4250000.0, 0.153), + ); + let parameters = Arc::new(CubicParameters::new_pure(propane).unwrap()); + let eos = Cubic::peng_robinson(parameters, None, None).unwrap(); + dbg!(&eos.critical_parameters); + let state = StateHD::new(300.0, 1e5, arr1(&[5.0])); + + let propane = PureRecord::new( + Identifier::new(None, Some("propane"), None, None, None, None), + 44.0962, + PengRobinsonRecord::new(369.96, 4250000.0, 0.153), + ); + let parameters = PengRobinsonParameters::new_pure(propane).unwrap(); + let pr = Arc::new(PengRobinson::new(Arc::new(parameters))); + + assert_eq!( + pr.residual_helmholtz_energy(&state), + eos.residual_helmholtz_energy(&state) + ) + } +} diff --git a/src/cubic/parameters.rs b/src/cubic/parameters.rs new file mode 100644 index 000000000..05301b9fd --- /dev/null +++ b/src/cubic/parameters.rs @@ -0,0 +1,217 @@ +use feos_core::parameter::{Identifier, Parameter, ParameterError, PureRecord}; +use ndarray::{Array1, Array2}; +use num_traits::Zero; +use serde::{Deserialize, Serialize}; + +/// Cubic parameters for a single substance. +#[derive(Serialize, Deserialize, Debug, Clone, Default)] +pub struct CubicRecord { + /// critical temperature in Kelvin + pub(crate) tc: f64, + /// critical pressure in Pascal + pub(crate) pc: f64, + /// acentric factor + pub(crate) acentric_factor: f64, +} + +impl CubicRecord { + /// Create a new pure substance record for the Cubic equation of state. + pub fn new(tc: f64, pc: f64, acentric_factor: f64) -> Self { + Self { + tc, + pc, + acentric_factor, + } + } +} + +impl std::fmt::Display for CubicRecord { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + write!(f, "CubicRecord(tc={} K", self.tc)?; + write!(f, ", pc={} Pa", self.pc)?; + write!(f, ", acentric factor={}", self.acentric_factor) + } +} + +/// Cubic binary interaction parameters. +#[derive(Serialize, Deserialize, Clone, Default)] +pub struct CubicBinaryRecord { + /// Binary interaction parameter for a + #[serde(skip_serializing_if = "f64::is_zero")] + #[serde(default)] + pub k_ij: f64, + /// Binary interaction parameter for b + #[serde(skip_serializing_if = "f64::is_zero")] + #[serde(default)] + pub l_ij: f64, + // /// Binary association parameters + // #[serde(flatten)] + // association: Option, +} + +impl CubicBinaryRecord { + pub fn new( + k_ij: Option, + l_ij: Option, + // rc_ab: Option, + // epsilon_k_ab: Option, + ) -> Self { + let k_ij = k_ij.unwrap_or_default(); + let l_ij = l_ij.unwrap_or_default(); + // let association = if rc_ab.is_none() && epsilon_k_ab.is_none() { + // None + // } else { + // Some(BinaryAssociationRecord::new(rc_ab, epsilon_k_ab, None)) + // }; + Self { + k_ij, + l_ij, + // association, + } + } +} + +impl From for CubicBinaryRecord { + fn from(k_ij: f64) -> Self { + Self { + k_ij, + l_ij: f64::default(), + // association: None, + } + } +} + +impl From for f64 { + fn from(binary_record: CubicBinaryRecord) -> Self { + binary_record.k_ij + } +} + +impl std::fmt::Display for CubicBinaryRecord { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + let mut tokens = vec![]; + if !self.k_ij.is_zero() { + tokens.push(format!("k_ij={}", self.k_ij)); + } + if !self.l_ij.is_zero() { + tokens.push(format!("l_ij={}", self.l_ij)); + } + // if let Some(association) = self.association { + // if let Some(rc_ab) = association.rc_ab { + // tokens.push(format!("rc_ab={}", rc_ab)); + // } + // if let Some(epsilon_k_ab) = association.epsilon_k_ab { + // tokens.push(format!("epsilon_k_ab={}", epsilon_k_ab)); + // } + // } + write!(f, "CubicBinaryRecord({})", tokens.join(", ")) + } +} + +/// Cubic parameters for one ore more substances. +pub struct CubicParameters { + /// Critical temperature in Kelvin + pub(super) tc: Array1, + pub(super) pc: Array1, + pub(super) acentric_factor: Array1, + /// Binary interaction parameter for a + pub(super) k_ij: Array2, + /// Binary interaction parameter for b + pub(super) l_ij: Array2, + /// Molar weight in units of g/mol + pub(super) molarweight: Array1, + /// List of pure component records + pub(super) pure_records: Vec>, + /// List of binary records + pub binary_records: Option>, +} + +impl std::fmt::Display for CubicParameters { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + self.pure_records + .iter() + .try_for_each(|pr| writeln!(f, "{}", pr))?; + writeln!(f, "\nk_ij:\n{}", self.k_ij) + } +} + +impl CubicParameters { + /// Build a simple parameter set without binary interaction parameters. + pub fn new_simple( + tc: &[f64], + pc: &[f64], + acentric_factor: &[f64], + molarweight: &[f64], + ) -> Result { + if [pc.len(), acentric_factor.len(), molarweight.len()] + .iter() + .any(|&l| l != tc.len()) + { + return Err(ParameterError::IncompatibleParameters(String::from( + "each component has to have parameters.", + ))); + } + let records = (0..tc.len()) + .map(|i| { + let record = CubicRecord { + tc: tc[i], + pc: pc[i], + acentric_factor: acentric_factor[i], + }; + let id = Identifier::default(); + PureRecord::new(id, molarweight[i], record) + }) + .collect(); + CubicParameters::from_records(records, None) + } +} + +impl Parameter for CubicParameters { + type Pure = CubicRecord; + type Binary = CubicBinaryRecord; + + /// Creates parameters from pure component records. + fn from_records( + pure_records: Vec>, + binary_records: Option>, + ) -> Result { + let n = pure_records.len(); + + let mut tc = Array1::zeros(n); + let mut pc = Array1::zeros(n); + let mut acentric_factor = Array1::zeros(n); + let mut molarweight = Array1::zeros(n); + + for (i, record) in pure_records.iter().enumerate() { + molarweight[i] = record.molarweight; + let r = &record.model_record; + tc[i] = r.tc; + pc[i] = r.pc; + acentric_factor[i] = r.acentric_factor; + } + + let br = binary_records.as_ref(); + let k_ij = br.map_or_else(|| Array2::zeros([n; 2]), |br| br.mapv(|br| br.k_ij)); + let l_ij = br.map_or_else(|| Array2::zeros([n; 2]), |br| br.mapv(|br| br.l_ij)); + + Ok(Self { + tc, + pc, + acentric_factor, + k_ij, + l_ij, + molarweight, + pure_records, + binary_records, + }) + } + + fn records( + &self, + ) -> ( + &[PureRecord], + Option<&Array2>, + ) { + (&self.pure_records, self.binary_records.as_ref()) + } +} diff --git a/src/cubic/python.rs b/src/cubic/python.rs new file mode 100644 index 000000000..9766434b5 --- /dev/null +++ b/src/cubic/python.rs @@ -0,0 +1,293 @@ +use super::alpha::{ + Alpha, GeneralizedTwu, MathiasCopeman, PengRobinson1976, PengRobinson1978, PengRobinson2019, + RedlichKwong1972, RedlichKwong2019, Soave, Twu, +}; +use super::mixing_rules::{MixingRule, Quadratic}; +use super::parameters::*; +use feos_core::parameter::*; +use feos_core::python::parameter::*; +use feos_core::{impl_binary_record, impl_json_handling, impl_parameter, impl_pure_record}; +use numpy::{prelude::*, PyReadonlyArray1}; +use numpy::{PyArray2, PyReadonlyArray2}; +use pyo3::exceptions::PyTypeError; +use pyo3::prelude::*; +use std::convert::{TryFrom, TryInto}; +use std::sync::Arc; + +#[pyclass(name = "Alpha")] +#[derive(Clone)] +pub struct PyAlpha(pub Alpha); + +#[pymethods] +impl PyAlpha { + /// Generalized Soave alpha function $\alpha(T) = \left[1 + m (1 - \sqrt{T_r})\right]^2$. + /// + /// with + /// + /// $m = \sum_i=0^N m_i \omega^i$ + /// + /// Parameters + /// ---------- + /// m : numpy.ndarray[float] + /// Polynomial coefficients for polynomial of acentric factors. + /// The number of coefficients specifies the polynomial degree. + /// + /// Returns + /// ------- + /// Alpha + #[staticmethod] + fn soave(m: PyReadonlyArray1<'_, f64>) -> Self { + Self(Soave::new(m.as_array().to_owned()).into()) + } + + /// Soave alpha function parameterized in 1972 for Redlich-Kwong EoS. + /// + /// $\alpha(T) = \left[1 + m_\text{RK} (1 - \sqrt{T_r})\right]^2$ + /// using + /// $m_\text{RK} = 0.480 + 1.574 \omega - 0.176 \omega^2$ + #[staticmethod] + fn redlich_kwing1972() -> Self { + Self(RedlichKwong1972.into()) + } + + /// Soave alpha function parameterized in 1976 for Peng Robinson EoS. + /// + /// $\alpha(T) = \left[1 + m_\text{PR} (1 - \sqrt{T_r})\right]^2$ + /// using + /// $m_\text{PR} = 0.37464 + 1.54226 \omega - 0.26992 \omega^2$ + #[staticmethod] + fn peng_robinson1976() -> Self { + Self(PengRobinson1976.into()) + } + + /// Soave alpha function parameterized in 1978 for Peng Robinson EoS. + /// + /// Uses different polynomials depending on acentric factor. + #[staticmethod] + fn peng_robinson1978() -> Self { + Self(PengRobinson1978.into()) + } + + /// Soave alpha function parameterized in 2019 for Redlich Kwong EoS. + /// + /// $\alpha(T) = \left[1 + m_\text{RK} (1 - \sqrt{T_r})\right]^2$ + /// using + /// $m_\text{RK} = x + x \omega - x \omega^2$ + #[staticmethod] + fn redlich_kwong2019() -> Self { + Self(RedlichKwong2019.into()) + } + + /// Soave alpha function parameterized in 2019 for Peng Robinson EoS. + /// + /// $\alpha(T) = \left[1 + m_\text{PR} (1 - \sqrt{T_r})\right]^2$ + /// using + /// $m_\text{PR} = x + x \omega - x \omega^2$ + #[staticmethod] + fn peng_robinson2019() -> Self { + Self(PengRobinson2019.into()) + } + + /// Mathias and Compeman alpha function (1983). + /// + /// $\alpha(T) = \left[1 + \sum_{i=1}^3 c_i (1 - \sqrt{T_r})^i\right]^2$ + /// with $c_i$ as substance specific adjustable parameters. + /// + /// Parameters + /// ---------- + /// c : List[[float; 3]] + /// The three `c` parameters for each substance. + /// + /// Returns + /// ------- + /// Alpha + #[staticmethod] + fn mathias_copeman(c: Vec<[f64; 3]>) -> Self { + Self(MathiasCopeman(c).into()) + } + + /// Generalized Twu alpha function (1995) parameterized for PR. + /// + /// Returns + /// ------- + /// Alpha + #[staticmethod] + fn peng_robinson_generalized_twu() -> Self { + Self(GeneralizedTwu::peng_robinson().into()) + } + + /// Generalized Twu alpha function (1995) parameterized for RK. + /// + /// Returns + /// ------- + /// Alpha + #[staticmethod] + fn redlich_kwong_generalized_twu() -> Self { + Self(GeneralizedTwu::redlich_kwong().into()) + } + + /// Generalized Twu alpha function (1991) parameterized for RK. + /// + /// If no `n` parameters are provided, the model of Twu (1988) is constructed. + /// + /// Parameters + /// ---------- + /// m : List[float] + /// l : List[float] + /// n : List[float], optional + /// Defaults to 2.0 for each component (Twu 1988). + /// + /// Returns + /// ------- + /// Alpha + #[staticmethod] + fn twu(l: Vec, m: Vec, n: Option>) -> Self { + Self(Twu::new(l, m, n).into()) + } +} + +/// Create a mixing rule +#[pyclass(name = "MixingRule")] +#[derive(Clone)] +pub struct PyMixingRule(pub MixingRule); + +#[pymethods] +impl PyMixingRule { + #[staticmethod] + fn quadratic() -> Self { + Self(Quadratic.into()) + } +} + +/// Create a set of cubic parameters from records. +#[pyclass(name = "CubicRecord")] +#[derive(Clone)] +pub struct PyCubicRecord(CubicRecord); + +#[pymethods] +impl PyCubicRecord { + #[new] + #[pyo3(text_signature = "(tc, pc, acentric_factor)")] + fn new(tc: f64, pc: f64, acentric_factor: f64) -> Self { + Self(CubicRecord::new(tc, pc, acentric_factor)) + } + + #[getter] + fn get_tc(&self) -> f64 { + self.0.tc + } + + #[getter] + fn get_pc(&self) -> f64 { + self.0.pc + } + + #[getter] + fn get_acentric_factor(&self) -> f64 { + self.0.acentric_factor + } + + fn __repr__(&self) -> PyResult { + Ok(self.0.to_string()) + } +} + +impl_json_handling!(PyCubicRecord); +impl_pure_record!(CubicRecord, PyCubicRecord); + +#[pyclass(name = "CubicBinaryRecord")] +#[derive(Clone)] +pub struct PyCubicBinaryRecord(CubicBinaryRecord); +impl_binary_record!(CubicBinaryRecord, PyCubicBinaryRecord); + +#[pyclass(name = "CubicParameters")] +#[derive(Clone)] +pub struct PyCubicParameters(pub Arc); + +#[pymethods] +impl PyCubicParameters { + // Create a set of cubic parameters from values. + /// + /// Parameters + /// ---------- + /// tc : float + /// critical temperature in units of Kelvin. + /// pc : float + /// critical pressure in units of Pascal. + /// acentric_factor: float + /// acentric factor + /// molarweight: float, optional + /// molar weight in units of Gram per Mol. + /// Returns + /// ------- + /// CubicParameters + #[pyo3(text_signature = "(tc, pc, acentric_factor, molarweight=None)")] + #[staticmethod] + fn from_values( + tc: f64, + pc: f64, + acentric_factor: f64, + molarweight: Option, + ) -> PyResult { + let pure_record = PureRecord::new( + Identifier::new( + Some(format!("{}", 1).as_str()), + None, + None, + None, + None, + None, + ), + molarweight.map_or(1.0, |v| v), + CubicRecord::new(tc, pc, acentric_factor), + ); + Ok(Self(Arc::new(CubicParameters::new_pure(pure_record)?))) + } + + #[getter] + fn get_k_ij<'py>(&self, py: Python<'py>) -> Option>> { + self.0 + .binary_records + .as_ref() + .map(|br| br.map(|br| br.k_ij).view().to_pyarray_bound(py)) + } + + #[getter] + fn get_l_ij<'py>(&self, py: Python<'py>) -> Option>> { + self.0 + .binary_records + .as_ref() + .map(|br| br.map(|br| br.l_ij).view().to_pyarray_bound(py)) + } + + fn _repr_markdown_(&self) -> String { + // self.0.to_markdown() + todo!() + } + + fn __repr__(&self) -> PyResult { + Ok(self.0.to_string()) + } +} + +impl_parameter!( + CubicParameters, + PyCubicParameters, + PyCubicRecord, + PyCubicBinaryRecord +); + +#[pymodule] +pub fn cubic(m: &Bound<'_, PyModule>) -> PyResult<()> { + m.add_class::()?; + m.add_class::()?; + m.add_class::()?; + + m.add_class::()?; + m.add_class::()?; + m.add_class::()?; + m.add_class::()?; + m.add_class::()?; + m.add_class::()?; + Ok(()) +} diff --git a/src/eos.rs b/src/eos.rs index 60bcf30fd..a659daa57 100644 --- a/src/eos.rs +++ b/src/eos.rs @@ -1,3 +1,5 @@ +#[cfg(feature = "cubic")] +use crate::cubic::Cubic; #[cfg(feature = "epcsaft")] use crate::epcsaft::ElectrolytePcSaft; #[cfg(feature = "gc_pcsaft")] @@ -48,4 +50,6 @@ pub enum ResidualModel { Pets(Pets), #[cfg(feature = "uvtheory")] UVTheory(UVTheory), + #[cfg(feature = "cubic")] + Cubic(Cubic), } diff --git a/src/python/cubic.rs b/src/python/cubic.rs index 93c7ee2af..8362de2fd 100644 --- a/src/python/cubic.rs +++ b/src/python/cubic.rs @@ -3,7 +3,7 @@ use feos_core::python::parameter::{PyChemicalRecord, PyIdentifier}; use pyo3::prelude::*; #[pymodule] -pub fn cubic(m: &Bound<'_, PyModule>) -> PyResult<()> { +pub fn cubic_old(m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_class::()?; m.add_class::()?; diff --git a/src/python/eos.rs b/src/python/eos.rs index 3494855c1..69b238087 100644 --- a/src/python/eos.rs +++ b/src/python/eos.rs @@ -1,3 +1,9 @@ +#[cfg(feature = "cubic")] +use crate::cubic::python::PyCubicParameters; +use crate::cubic::python::{PyAlpha, PyMixingRule}; +#[cfg(feature = "cubic")] +use crate::cubic::Cubic; +use crate::cubic::CubicOptions; use crate::eos::ResidualModel; #[cfg(feature = "epcsaft")] use crate::epcsaft::python::PyElectrolytePcSaftParameters; @@ -240,6 +246,89 @@ impl PyEquationOfState { Self(Arc::new(EquationOfState::new(ideal_gas, residual))) } + /// Generalized cubic equations of state. + /// + /// Parameters + /// ---------- + /// parameters : CubicParameters + /// The parameters of the equation of state to use. + /// + /// Returns + /// ------- + /// EquationOfState + /// The equation of state that can be used to compute thermodynamic + /// states. + #[staticmethod] + pub fn cubic( + parameters: PyCubicParameters, + alpha: PyAlpha, + mixing_rule: PyMixingRule, + delta: (f64, f64), + ) -> PyResult { + let options = CubicOptions { + alpha: alpha.0, + mixing: mixing_rule.0, + delta: delta.into(), + }; + let residual = Arc::new(ResidualModel::Cubic(Cubic::new(parameters.0, options)?)); + let ideal_gas = Arc::new(IdealGasModel::NoModel(residual.components())); + Ok(Self(Arc::new(EquationOfState::new(ideal_gas, residual)))) + } + + /// Peng-Robinson equations of state. + /// + /// Parameters + /// ---------- + /// parameters : CubicParameters + /// The parameters of the equation of state to use. + /// + /// Returns + /// ------- + /// EquationOfState + /// The equation of state that can be used to compute thermodynamic + /// states. + #[staticmethod] + pub fn peng_robinson( + parameters: PyCubicParameters, + alpha: Option, + mixing_rule: Option, + ) -> PyResult { + let residual = Arc::new(ResidualModel::Cubic(Cubic::peng_robinson( + parameters.0, + alpha.map(|a| a.0), + mixing_rule.map(|mr| mr.0), + )?)); + let ideal_gas = Arc::new(IdealGasModel::NoModel(residual.components())); + Ok(Self(Arc::new(EquationOfState::new(ideal_gas, residual)))) + } + + /// Redlich-Kwong equations of state. + /// + /// Parameters + /// ---------- + /// parameters : CubicParameters + /// The parameters of the equation of state to use. + /// + /// Returns + /// ------- + /// EquationOfState + /// The equation of state that can be used to compute thermodynamic + /// states. + #[staticmethod] + pub fn redlich_kwong( + parameters: PyCubicParameters, + alpha: Option, + mixing_rule: Option, + ) -> PyResult { + let residual = Arc::new(ResidualModel::Cubic(Cubic::redlich_kwong( + parameters.0, + alpha.map(|a| a.0), + mixing_rule.map(|mr| mr.0), + )?)); + let ideal_gas = Arc::new(IdealGasModel::NoModel(residual.components())); + Ok(Self(Arc::new(EquationOfState::new(ideal_gas, residual)))) + } + /// Peng-Robinson equation of state. /// /// Parameters @@ -253,7 +342,7 @@ impl PyEquationOfState { /// The PR equation of state that can be used to compute thermodynamic /// states. #[staticmethod] - pub fn peng_robinson(parameters: PyPengRobinsonParameters) -> Self { + pub fn peng_robinson_old(parameters: PyPengRobinsonParameters) -> Self { let residual = Arc::new(ResidualModel::PengRobinson(PengRobinson::new(parameters.0))); let ideal_gas = Arc::new(IdealGasModel::NoModel(residual.components())); Self(Arc::new(EquationOfState::new(ideal_gas, residual))) diff --git a/src/python/mod.rs b/src/python/mod.rs index e1613b673..f4a10f51e 100644 --- a/src/python/mod.rs +++ b/src/python/mod.rs @@ -1,3 +1,5 @@ +#[cfg(feature = "cubic")] +use crate::cubic::python::cubic as cubic_module; #[cfg(feature = "epcsaft")] use crate::epcsaft::python::epcsaft as epcsaft_module; #[cfg(feature = "gc_pcsaft")] @@ -21,7 +23,7 @@ mod cubic; mod dippr; mod eos; mod joback; -use cubic::cubic as cubic_module; +use cubic::cubic_old as cubic_module_old; use dippr::dippr as dippr_module; use eos::eos as eos_module; use joback::joback as joback_module; @@ -41,6 +43,8 @@ pub fn feos(m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_wrapped(wrap_pymodule!(dft_module))?; m.add_wrapped(wrap_pymodule!(joback_module))?; m.add_wrapped(wrap_pymodule!(dippr_module))?; + m.add_wrapped(wrap_pymodule!(cubic_module_old))?; + #[cfg(feature = "cubic")] m.add_wrapped(wrap_pymodule!(cubic_module))?; #[cfg(feature = "pcsaft")] m.add_wrapped(wrap_pymodule!(pcsaft_module))?; @@ -67,6 +71,7 @@ pub fn feos(m: &Bound<'_, PyModule>) -> PyResult<()> { set_path(m, "feos.dft.estimator", "dft.estimator_dft")?; set_path(m, "feos.joback", "joback")?; set_path(m, "feos.dippr", "dippr")?; + set_path(m, "feos.cubic_old", "cubic_old")?; set_path(m, "feos.cubic", "cubic")?; #[cfg(feature = "pcsaft")] set_path(m, "feos.pcsaft", "pcsaft")?; From e44d857ffdf09068a7a39106cebd796abadb3afa Mon Sep 17 00:00:00 2001 From: Gernot Bauer Date: Fri, 7 Jun 2024 15:58:51 +0200 Subject: [PATCH 3/4] fixed typo --- feos-core/src/cubic.rs | 2 -- src/cubic/python.rs | 2 +- 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/feos-core/src/cubic.rs b/feos-core/src/cubic.rs index bb51100a3..1b178ef51 100644 --- a/feos-core/src/cubic.rs +++ b/feos-core/src/cubic.rs @@ -200,8 +200,6 @@ impl Residual for PengRobinson { // Helmholtz energy let n = state.moles.sum(); let v = state.volume; - dbg!(&b); - dbg!(&ak_mix); n * ((v / (v - b * n)).ln() - ak_mix / (b * SQRT_2 * 2.0 * state.temperature) * ((v + b * n * (1.0 + SQRT_2)) / (v + b * n * (1.0 - SQRT_2))).ln()) diff --git a/src/cubic/python.rs b/src/cubic/python.rs index 9766434b5..9f18b6f06 100644 --- a/src/cubic/python.rs +++ b/src/cubic/python.rs @@ -46,7 +46,7 @@ impl PyAlpha { /// using /// $m_\text{RK} = 0.480 + 1.574 \omega - 0.176 \omega^2$ #[staticmethod] - fn redlich_kwing1972() -> Self { + fn redlich_kwong1972() -> Self { Self(RedlichKwong1972.into()) } From 9b5df640be43ee438123baef5265bfd79504df89 Mon Sep 17 00:00:00 2001 From: Gernot Bauer Date: Fri, 7 Jun 2024 16:08:08 +0200 Subject: [PATCH 4/4] fixed bug in subset method --- src/cubic/mod.rs | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/src/cubic/mod.rs b/src/cubic/mod.rs index 6f33bdc63..a696d8c64 100644 --- a/src/cubic/mod.rs +++ b/src/cubic/mod.rs @@ -104,6 +104,16 @@ pub struct CubicOptions { pub(crate) delta: Delta, } +impl CubicOptions { + fn subset(&self, component_list: &[usize]) -> Self { + Self { + alpha: self.alpha.subset(component_list), + mixing: self.mixing.clone(), + delta: self.delta.clone(), + } + } +} + /// A generic cubic equation of state. pub struct Cubic { /// Parameters @@ -198,7 +208,7 @@ impl Components for Cubic { fn subset(&self, component_list: &[usize]) -> Self { Self { parameters: Arc::new(self.parameters.subset(component_list)), - options: self.options.clone(), + options: self.options.subset(component_list), critical_parameters: self.critical_parameters.subset(component_list), } }