diff --git a/feos-dft/CHANGELOG.md b/feos-dft/CHANGELOG.md index 8cf352dc0..387ed5e4f 100644 --- a/feos-dft/CHANGELOG.md +++ b/feos-dft/CHANGELOG.md @@ -5,6 +5,9 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). ## [Unreleased] +### Added +- `PlanarInterface` now has methods for the calculation of the 90-10 interface thickness (`PlanarInterface::interfacial_thickness`), interfacial enrichtment (`PlanarInterface::interfacial_enrichment`), and relative adsorption (`PlanarInterface::relative_adsorption`). + ### Changed - Added `Send` and `Sync` as supertraits to `HelmholtzEnergyFunctional` and all related traits. [#57](https://github.com/feos-org/feos/pull/57) - Renamed the `3d_dft` feature to `rayon` in accordance to the other feos crates. [#62](https://github.com/feos-org/feos/pull/62) diff --git a/feos-dft/src/interface/mod.rs b/feos-dft/src/interface/mod.rs index 0c9841790..956288327 100644 --- a/feos-dft/src/interface/mod.rs +++ b/feos-dft/src/interface/mod.rs @@ -196,6 +196,140 @@ impl PlanarInterface { self } + /// Relative adsorption of component `i' with respect to `j': \Gamma_i^(j) + pub fn relative_adsorption(&self) -> EosResult> { + let s = self.profile.density.shape(); + let mut rho_l = Array1::zeros(s[0]) * U::reference_density(); + let mut rho_v = Array1::zeros(s[0]) * U::reference_density(); + + // Calculate the partial densities in the liquid and in the vapor phase + for i in 0..s[0] { + rho_l.try_set(i, self.profile.density.get((i, 0)))?; + rho_v.try_set(i, self.profile.density.get((i, s[1] - 1)))?; + } + + // Calculate \Gamma_i^(j) + let gamma = QuantityArray2::from_shape_fn((s[0], s[0]), |(i, j)| { + if i == j { + 0.0 * U::reference_density() * U::reference_length() + } else { + -(rho_l.get(i) - rho_v.get(i)) + * ((&self.profile.density.index_axis(Axis_nd(0), j) - rho_l.get(j)) + / (rho_l.get(j) - rho_v.get(j)) + - (&self.profile.density.index_axis(Axis_nd(0), i) - rho_l.get(i)) + / (rho_l.get(i) - rho_v.get(i))) + .integrate(&self.profile.grid.integration_weights_unit()) + } + }); + + Ok(gamma) + } + + /// Interfacial enrichment of component `i': E_i + pub fn interfacial_enrichment(&self) -> EosResult> { + let s = self.profile.density.shape(); + let density = self.profile.density.to_reduced(U::reference_density())?; + let rho_l = density.index_axis(Axis_nd(1), 0); + let rho_v = density.index_axis(Axis_nd(1), s[1] - 1); + + let enrichment = Array1::from_shape_fn(s[0], |i| { + *(density + .index_axis(Axis_nd(0), i) + .iter() + .max_by(|&a, &b| a.total_cmp(b)) + .unwrap()) // panics only of iterator is empty + / rho_l[i].max(rho_v[i]) + }); + + Ok(enrichment) + } + + /// Interface thickness (90-10 number density difference) + pub fn interfacial_thickness(&self) -> EosResult> { + let s = self.profile.density.shape(); + let rho = self + .profile + .density + .sum_axis(Axis_nd(0)) + .to_reduced(U::reference_density())?; + let z = self.profile.grid.grids()[0]; + let dz = z[1] - z[0]; + + let limits = (0.9_f64, 0.1_f64); + let (limit_upper, limit_lower) = if limits.0 > limits.1 { + (limits.0, limits.1) + } else { + (limits.1, limits.0) + }; + + if limit_upper >= 1.0 || limit_upper.is_sign_negative() { + return Err(EosError::IterationFailed(String::from( + "Upper limit 'l' of interface thickness needs to satisfy 0 < l < 1.", + ))); + } + if limit_lower >= 1.0 || limit_lower.is_sign_negative() { + return Err(EosError::IterationFailed(String::from( + "Lower limit 'l' of interface thickness needs to satisfy 0 < l < 1.", + ))); + } + + // Get the densities in the liquid and in the vapor phase + let rho_v = rho[0].min(rho[s[1] - 1]); + let rho_l = rho[0].max(rho[s[1] - 1]); + + if (rho_l - rho_v).abs() < 1.0e-10 { + return Ok(0.0 * U::reference_length()); + } + + // Density boundaries for interface definition + let rho_upper = rho_v + limit_upper * (rho_l - rho_v); + let rho_lower = rho_v + limit_lower * (rho_l - rho_v); + + // Get indizes right of intersection between density profile and + // constant density boundaries + let index_upper_plus = if rho[0] >= rho[s[1] - 1] { + rho.iter() + .enumerate() + .find(|(_, &x)| (x - rho_upper).is_sign_negative()) + .expect("Could not find rho_upper value!") + .0 + } else { + rho.iter() + .enumerate() + .find(|(_, &x)| (rho_upper - x).is_sign_negative()) + .expect("Could not find rho_upper value!") + .0 + }; + let index_lower_plus = if rho[0] >= rho[s[1] - 1] { + rho.iter() + .enumerate() + .find(|(_, &x)| (x - rho_lower).is_sign_negative()) + .expect("Could not find rho_lower value!") + .0 + } else { + rho.iter() + .enumerate() + .find(|(_, &x)| (rho_lower - x).is_sign_negative()) + .expect("Could not find rho_lower value!") + .0 + }; + + // Calculate distance between two density points using a linear + // interpolated density profiles between the two grid points where the + // density profile crosses the limiting densities + let z_upper = z[index_upper_plus - 1] + + (rho_upper - rho[index_upper_plus - 1]) + / (rho[index_upper_plus] - rho[index_upper_plus - 1]) + * dz; + let z_lower = z[index_lower_plus - 1] + + (rho_lower - rho[index_lower_plus - 1]) + / (rho[index_lower_plus] - rho[index_lower_plus - 1]) + * dz; + + // Return + Ok((z_lower - z_upper) * U::reference_length()) + } + fn set_density_scale(&mut self, init: &QuantityArray2) { assert_eq!(self.profile.density.shape(), init.shape()); let n_grid = self.profile.density.shape()[1]; diff --git a/feos-dft/src/interface/surface_tension_diagram.rs b/feos-dft/src/interface/surface_tension_diagram.rs index 0122432ee..84b768fc1 100644 --- a/feos-dft/src/interface/surface_tension_diagram.rs +++ b/feos-dft/src/interface/surface_tension_diagram.rs @@ -2,7 +2,8 @@ use super::PlanarInterface; use crate::functional::{HelmholtzEnergyFunctional, DFT}; use crate::solver::DFTSolver; use feos_core::{EosUnit, PhaseEquilibrium, StateVec}; -use quantity::{QuantityArray1, QuantityScalar}; +use ndarray::Array1; +use quantity::{QuantityArray1, QuantityArray2, QuantityScalar}; const DEFAULT_GRID_POINTS: usize = 2048; @@ -76,4 +77,25 @@ impl SurfaceTensionDiagram { self.profiles[i].surface_tension.unwrap() }) } + + pub fn relative_adsorption(&self) -> Vec> { + self.profiles + .iter() + .map(|planar_interf| planar_interf.relative_adsorption().unwrap()) + .collect() + } + + pub fn interfacial_enrichment(&self) -> Vec> { + self.profiles + .iter() + .map(|planar_interf| planar_interf.interfacial_enrichment().unwrap()) + .collect() + } + + pub fn interfacial_thickness(&self) -> QuantityArray1 { + self.profiles + .iter() + .map(|planar_interf| planar_interf.interfacial_thickness().unwrap()) + .collect() + } } diff --git a/feos-dft/src/python/interface/mod.rs b/feos-dft/src/python/interface/mod.rs index f0787b4ab..fc4ae5df8 100644 --- a/feos-dft/src/python/interface/mod.rs +++ b/feos-dft/src/python/interface/mod.rs @@ -114,5 +114,42 @@ macro_rules! impl_planar_interface { PyPhaseEquilibrium(self.0.vle.clone()) } } + + #[pymethods] + impl PyPlanarInterface { + /// Calculates the Gibbs' relative adsorption of component `i' with + /// respect to `j': \Gamma_i^(j) + /// + /// Returns + /// ------- + /// SIArray2 + /// + #[pyo3(text_signature = "($self)")] + fn relative_adsorption(&self) -> PyResult { + Ok(self.0.relative_adsorption()?.into()) + } + + /// Calculates the interfacial enrichment E_i. + /// + /// Returns + /// ------- + /// numpy.ndarray + /// + #[pyo3(text_signature = "($self)")] + fn interfacial_enrichment<'py>(&self, py: Python<'py>) -> PyResult<&'py PyArray1> { + Ok(self.0.interfacial_enrichment()?.to_pyarray(py)) + } + + /// Calculates the interfacial thickness (90-10 number density difference) + /// + /// Returns + /// ------- + /// SINumber + /// + #[pyo3(text_signature = "($self)")] + fn interfacial_thickness(&self) -> PyResult { + Ok(self.0.interfacial_thickness()?.into()) + } + } }; } diff --git a/feos-dft/src/python/interface/surface_tension_diagram.rs b/feos-dft/src/python/interface/surface_tension_diagram.rs index 0f2687b8a..eaaafc602 100644 --- a/feos-dft/src/python/interface/surface_tension_diagram.rs +++ b/feos-dft/src/python/interface/surface_tension_diagram.rs @@ -75,6 +75,21 @@ macro_rules! impl_surface_tension_diagram { pub fn get_surface_tension(&mut self) -> PySIArray1 { self.0.surface_tension().into() } + + #[getter] + pub fn get_relative_adsorption(&self) -> Vec { + self.0.relative_adsorption().iter().cloned().map(|i| i.into()).collect() + } + + #[getter] + pub fn get_interfacial_enrichment<'py>(&self, py: Python<'py>) -> Vec<&'py PyArray1> { + self.0.interfacial_enrichment().iter().map(|i| i.to_pyarray(py)).collect() + } + + #[getter] + pub fn interfacial_thickness(&self) -> PySIArray1 { + self.0.interfacial_thickness().into() + } } }; }