Skip to content

Commit 8f0ebcc

Browse files
RolfStierleRolf Stierle
RolfStierle
and
Rolf Stierle
authored
Relative adsorption, interfacial enrichment, and interface thickness. (#64)
* Relative adsorption, interfacial enrichment, and interface thickness (90-10; just in Python interface, otherwise general) implemented. * Improvements for Interfacial property diagrams. Co-authored-by: Rolf Stierle <stierle@itt.uni-stuttgart.de>
1 parent 329770b commit 8f0ebcc

File tree

5 files changed

+212
-1
lines changed

5 files changed

+212
-1
lines changed

‎feos-dft/CHANGELOG.md‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,9 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
55
and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
66

77
## [Unreleased]
8+
### Added
9+
- `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`).
10+
811
### Changed
912
- Added `Send` and `Sync` as supertraits to `HelmholtzEnergyFunctional` and all related traits. [#57](https://github.com/feos-org/feos/pull/57)
1013
- Renamed the `3d_dft` feature to `rayon` in accordance to the other feos crates. [#62](https://github.com/feos-org/feos/pull/62)

‎feos-dft/src/interface/mod.rs‎

Lines changed: 134 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -205,6 +205,140 @@ impl<U: EosUnit, F: HelmholtzEnergyFunctional> PlanarInterface<U, F> {
205205
self
206206
}
207207

208+
/// Relative adsorption of component `i' with respect to `j': \Gamma_i^(j)
209+
pub fn relative_adsorption(&self) -> EosResult<QuantityArray2<U>> {
210+
let s = self.profile.density.shape();
211+
let mut rho_l = Array1::zeros(s[0]) * U::reference_density();
212+
let mut rho_v = Array1::zeros(s[0]) * U::reference_density();
213+
214+
// Calculate the partial densities in the liquid and in the vapor phase
215+
for i in 0..s[0] {
216+
rho_l.try_set(i, self.profile.density.get((i, 0)))?;
217+
rho_v.try_set(i, self.profile.density.get((i, s[1] - 1)))?;
218+
}
219+
220+
// Calculate \Gamma_i^(j)
221+
let gamma = QuantityArray2::from_shape_fn((s[0], s[0]), |(i, j)| {
222+
if i == j {
223+
0.0 * U::reference_density() * U::reference_length()
224+
} else {
225+
-(rho_l.get(i) - rho_v.get(i))
226+
* ((&self.profile.density.index_axis(Axis_nd(0), j) - rho_l.get(j))
227+
/ (rho_l.get(j) - rho_v.get(j))
228+
- (&self.profile.density.index_axis(Axis_nd(0), i) - rho_l.get(i))
229+
/ (rho_l.get(i) - rho_v.get(i)))
230+
.integrate(&self.profile.grid.integration_weights_unit())
231+
}
232+
});
233+
234+
Ok(gamma)
235+
}
236+
237+
/// Interfacial enrichment of component `i': E_i
238+
pub fn interfacial_enrichment(&self) -> EosResult<Array1<f64>> {
239+
let s = self.profile.density.shape();
240+
let density = self.profile.density.to_reduced(U::reference_density())?;
241+
let rho_l = density.index_axis(Axis_nd(1), 0);
242+
let rho_v = density.index_axis(Axis_nd(1), s[1] - 1);
243+
244+
let enrichment = Array1::from_shape_fn(s[0], |i| {
245+
*(density
246+
.index_axis(Axis_nd(0), i)
247+
.iter()
248+
.max_by(|&a, &b| a.total_cmp(b))
249+
.unwrap()) // panics only of iterator is empty
250+
/ rho_l[i].max(rho_v[i])
251+
});
252+
253+
Ok(enrichment)
254+
}
255+
256+
/// Interface thickness (90-10 number density difference)
257+
pub fn interfacial_thickness(&self) -> EosResult<QuantityScalar<U>> {
258+
let s = self.profile.density.shape();
259+
let rho = self
260+
.profile
261+
.density
262+
.sum_axis(Axis_nd(0))
263+
.to_reduced(U::reference_density())?;
264+
let z = self.profile.grid.grids()[0];
265+
let dz = z[1] - z[0];
266+
267+
let limits = (0.9_f64, 0.1_f64);
268+
let (limit_upper, limit_lower) = if limits.0 > limits.1 {
269+
(limits.0, limits.1)
270+
} else {
271+
(limits.1, limits.0)
272+
};
273+
274+
if limit_upper >= 1.0 || limit_upper.is_sign_negative() {
275+
return Err(EosError::IterationFailed(String::from(
276+
"Upper limit 'l' of interface thickness needs to satisfy 0 < l < 1.",
277+
)));
278+
}
279+
if limit_lower >= 1.0 || limit_lower.is_sign_negative() {
280+
return Err(EosError::IterationFailed(String::from(
281+
"Lower limit 'l' of interface thickness needs to satisfy 0 < l < 1.",
282+
)));
283+
}
284+
285+
// Get the densities in the liquid and in the vapor phase
286+
let rho_v = rho[0].min(rho[s[1] - 1]);
287+
let rho_l = rho[0].max(rho[s[1] - 1]);
288+
289+
if (rho_l - rho_v).abs() < 1.0e-10 {
290+
return Ok(0.0 * U::reference_length());
291+
}
292+
293+
// Density boundaries for interface definition
294+
let rho_upper = rho_v + limit_upper * (rho_l - rho_v);
295+
let rho_lower = rho_v + limit_lower * (rho_l - rho_v);
296+
297+
// Get indizes right of intersection between density profile and
298+
// constant density boundaries
299+
let index_upper_plus = if rho[0] >= rho[s[1] - 1] {
300+
rho.iter()
301+
.enumerate()
302+
.find(|(_, &x)| (x - rho_upper).is_sign_negative())
303+
.expect("Could not find rho_upper value!")
304+
.0
305+
} else {
306+
rho.iter()
307+
.enumerate()
308+
.find(|(_, &x)| (rho_upper - x).is_sign_negative())
309+
.expect("Could not find rho_upper value!")
310+
.0
311+
};
312+
let index_lower_plus = if rho[0] >= rho[s[1] - 1] {
313+
rho.iter()
314+
.enumerate()
315+
.find(|(_, &x)| (x - rho_lower).is_sign_negative())
316+
.expect("Could not find rho_lower value!")
317+
.0
318+
} else {
319+
rho.iter()
320+
.enumerate()
321+
.find(|(_, &x)| (rho_lower - x).is_sign_negative())
322+
.expect("Could not find rho_lower value!")
323+
.0
324+
};
325+
326+
// Calculate distance between two density points using a linear
327+
// interpolated density profiles between the two grid points where the
328+
// density profile crosses the limiting densities
329+
let z_upper = z[index_upper_plus - 1]
330+
+ (rho_upper - rho[index_upper_plus - 1])
331+
/ (rho[index_upper_plus] - rho[index_upper_plus - 1])
332+
* dz;
333+
let z_lower = z[index_lower_plus - 1]
334+
+ (rho_lower - rho[index_lower_plus - 1])
335+
/ (rho[index_lower_plus] - rho[index_lower_plus - 1])
336+
* dz;
337+
338+
// Return
339+
Ok((z_lower - z_upper) * U::reference_length())
340+
}
341+
208342
fn set_density_scale(&mut self, init: &QuantityArray2<U>) {
209343
assert_eq!(self.profile.density.shape(), init.shape());
210344
let n_grid = self.profile.density.shape()[1];

‎feos-dft/src/interface/surface_tension_diagram.rs‎

Lines changed: 23 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,8 @@ use super::PlanarInterface;
22
use crate::functional::{HelmholtzEnergyFunctional, DFT};
33
use crate::solver::DFTSolver;
44
use feos_core::{EosUnit, PhaseEquilibrium, StateVec};
5-
use quantity::{QuantityArray1, QuantityScalar};
5+
use ndarray::Array1;
6+
use quantity::{QuantityArray1, QuantityArray2, QuantityScalar};
67

78
const DEFAULT_GRID_POINTS: usize = 2048;
89

@@ -79,4 +80,25 @@ impl<U: EosUnit, F: HelmholtzEnergyFunctional> SurfaceTensionDiagram<U, F> {
7980
self.profiles[i].surface_tension.unwrap()
8081
})
8182
}
83+
84+
pub fn relative_adsorption(&self) -> Vec<QuantityArray2<U>> {
85+
self.profiles
86+
.iter()
87+
.map(|planar_interf| planar_interf.relative_adsorption().unwrap())
88+
.collect()
89+
}
90+
91+
pub fn interfacial_enrichment(&self) -> Vec<Array1<f64>> {
92+
self.profiles
93+
.iter()
94+
.map(|planar_interf| planar_interf.interfacial_enrichment().unwrap())
95+
.collect()
96+
}
97+
98+
pub fn interfacial_thickness(&self) -> QuantityArray1<U> {
99+
self.profiles
100+
.iter()
101+
.map(|planar_interf| planar_interf.interfacial_thickness().unwrap())
102+
.collect()
103+
}
82104
}

‎feos-dft/src/python/interface/mod.rs‎

Lines changed: 37 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -124,5 +124,42 @@ macro_rules! impl_planar_interface {
124124
PyPhaseEquilibrium(self.0.vle.clone())
125125
}
126126
}
127+
128+
#[pymethods]
129+
impl PyPlanarInterface {
130+
/// Calculates the Gibbs' relative adsorption of component `i' with
131+
/// respect to `j': \Gamma_i^(j)
132+
///
133+
/// Returns
134+
/// -------
135+
/// SIArray2
136+
///
137+
#[pyo3(text_signature = "($self)")]
138+
fn relative_adsorption(&self) -> PyResult<PySIArray2> {
139+
Ok(self.0.relative_adsorption()?.into())
140+
}
141+
142+
/// Calculates the interfacial enrichment E_i.
143+
///
144+
/// Returns
145+
/// -------
146+
/// numpy.ndarray
147+
///
148+
#[pyo3(text_signature = "($self)")]
149+
fn interfacial_enrichment<'py>(&self, py: Python<'py>) -> PyResult<&'py PyArray1<f64>> {
150+
Ok(self.0.interfacial_enrichment()?.to_pyarray(py))
151+
}
152+
153+
/// Calculates the interfacial thickness (90-10 number density difference)
154+
///
155+
/// Returns
156+
/// -------
157+
/// SINumber
158+
///
159+
#[pyo3(text_signature = "($self)")]
160+
fn interfacial_thickness(&self) -> PyResult<PySINumber> {
161+
Ok(self.0.interfacial_thickness()?.into())
162+
}
163+
}
127164
};
128165
}

‎feos-dft/src/python/interface/surface_tension_diagram.rs‎

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -81,6 +81,21 @@ macro_rules! impl_surface_tension_diagram {
8181
pub fn get_surface_tension(&mut self) -> PySIArray1 {
8282
self.0.surface_tension().into()
8383
}
84+
85+
#[getter]
86+
pub fn get_relative_adsorption(&self) -> Vec<PySIArray2> {
87+
self.0.relative_adsorption().iter().cloned().map(|i| i.into()).collect()
88+
}
89+
90+
#[getter]
91+
pub fn get_interfacial_enrichment<'py>(&self, py: Python<'py>) -> Vec<&'py PyArray1<f64>> {
92+
self.0.interfacial_enrichment().iter().map(|i| i.to_pyarray(py)).collect()
93+
}
94+
95+
#[getter]
96+
pub fn interfacial_thickness(&self) -> PySIArray1 {
97+
self.0.interfacial_thickness().into()
98+
}
8499
}
85100
};
86101
}

0 commit comments

Comments
 (0)