Revive phi in hetero gc PC-SAFT (#157)

feos-org · prehner · Apr 14, 2023 · May 30, 2023 · May 31, 2023 · Jul 7, 2023
commit 16787018882e52e622ad33ebb1fdf7687f471cdb
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -8,6 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Changed
 - Changed the internal implementation of the association contribution to accomodate more general association schemes. [#150](https://github.com/feos-org/feos/pull/150)
 - To comply with the new association implementation, the default values of `na` and `nb` are now `0` rather than `1`. Parameter files have been adapted accordingly. [#150](https://github.com/feos-org/feos/pull/150)
+- Added the possibility to specify a pure component correction parameter `phi` for the heterosegmented gc PC-SAFT equation of state. [#157](https://github.com/feos-org/feos/pull/157)
 
 ## [0.4.3] - 2023-03-20
 - Python only: Release the changes introduced in `feos-core` 0.4.2.

diff --git a/src/gc_pcsaft/eos/parameter.rs b/src/gc_pcsaft/eos/parameter.rs
@@ -19,18 +19,21 @@ pub struct GcPcSaftChemicalRecord {
     pub identifier: Identifier,
     pub segments: HashMap<String, f64>,
     pub bonds: HashMap<[String; 2], f64>,
+    phi: f64,
 }
 
 impl GcPcSaftChemicalRecord {
     pub fn new(
         identifier: Identifier,
         segments: HashMap<String, f64>,
         bonds: HashMap<[String; 2], f64>,
+        phi: f64,
     ) -> Self {
         Self {
             identifier,
             segments,
             bonds,
+            phi,
         }
     }
 }
@@ -53,11 +56,13 @@ impl From<ChemicalRecord> for GcPcSaftChemicalRecord {
             chemical_record.identifier.clone(),
             chemical_record.segment_count(),
             chemical_record.bond_count(),
+            1.0,
         )
     }
 }
 
 /// Parameter set required for the gc-PC-SAFT equation of state.
+#[derive(Clone)]
 pub struct GcPcSaftEosParameters {
     pub molarweight: Array1<f64>,
     pub component_index: Array1<usize>,
@@ -116,11 +121,14 @@ impl ParameterHetero for GcPcSaftEosParameters {
         let mut sigma_mix = Vec::new();
         let mut epsilon_k_mix = Vec::new();
 
+        let mut phi = Vec::new();
+
         let mut joback_records = Vec::new();
 
         for (i, chemical_record) in chemical_records.iter().cloned().enumerate() {
             let mut segment_indices = IndexMap::with_capacity(segment_records.len());
             let segment_map = chemical_record.segment_map(&segment_records)?;
+            phi.push(chemical_record.phi);
 
             let mut m_i = 0.0;
             let mut sigma_i = 0.0;
@@ -214,7 +222,7 @@ impl ParameterHetero for GcPcSaftEosParameters {
         let sigma_ij =
             Array2::from_shape_fn([sigma.len(); 2], |(i, j)| 0.5 * (sigma[i] + sigma[j]));
         let epsilon_k_ij = Array2::from_shape_fn([epsilon_k.len(); 2], |(i, j)| {
-            (epsilon_k[i] * epsilon_k[j]).sqrt()
+            (epsilon_k[i] * phi[component_index[i]] * epsilon_k[j] * phi[component_index[j]]).sqrt()
         }) * (1.0 - &k_ij);
 
         // Combining rules polar
@@ -271,6 +279,14 @@ impl ParameterHetero for GcPcSaftEosParameters {
     }
 }
 
+impl GcPcSaftEosParameters {
+    pub fn phi(self, phi: &[f64]) -> Result<Self, ParameterError> {
+        let mut cr = self.chemical_records;
+        cr.iter_mut().zip(phi.iter()).for_each(|(c, &p)| c.phi = p);
+        Self::from_segments(cr, self.segment_records, self.binary_segment_records)
+    }
+}
+
 impl HardSphereProperties for GcPcSaftEosParameters {
     fn monomer_shape<N: DualNum<f64>>(&self, _: N) -> MonomerShape<N> {
         let m = self.m.mapv(N::from);

diff --git a/src/gc_pcsaft/python/mod.rs b/src/gc_pcsaft/python/mod.rs
@@ -95,6 +95,10 @@ impl_parameter_from_segments!(GcPcSaftEosParameters, PyGcPcSaftEosParameters);
 
 #[pymethods]
 impl PyGcPcSaftEosParameters {
+    fn phi(&self, phi: Vec<f64>) -> PyResult<Self> {
+        Ok(Self(Arc::new((*self.0).clone().phi(&phi)?)))
+    }
+
     fn _repr_markdown_(&self) -> String {
         self.0.to_markdown()
     }