feos_dft/interface/
mod.rs

1//! Density profiles at planar interfaces and interfacial tensions.
2use crate::functional::HelmholtzEnergyFunctional;
3use crate::geometry::{Axis, Grid};
4use crate::pdgt::PdgtFunctionalProperties;
5use crate::profile::{DFTProfile, DFTSpecifications};
6use crate::solver::DFTSolver;
7use feos_core::{Contributions, FeosError, FeosResult, PhaseEquilibrium, ReferenceSystem};
8use ndarray::{Array1, Array2, Axis as Axis_nd, Ix1, s};
9use quantity::{Area, Density, Length, Moles, SurfaceTension, Temperature};
10use std::sync::Arc;
11
12mod surface_tension_diagram;
13pub use surface_tension_diagram::SurfaceTensionDiagram;
14
15const RELATIVE_WIDTH: f64 = 6.0;
16const MIN_WIDTH: f64 = 100.0;
17
18/// Density profile and properties of a planar interface.
19#[derive(Clone)]
20pub struct PlanarInterface<F: HelmholtzEnergyFunctional> {
21    pub profile: DFTProfile<Ix1, F>,
22    pub vle: PhaseEquilibrium<F, 2>,
23    pub surface_tension: Option<SurfaceTension>,
24    pub equimolar_radius: Option<Length>,
25}
26
27impl<F: HelmholtzEnergyFunctional> PlanarInterface<F> {
28    pub fn solve_inplace(&mut self, solver: Option<&DFTSolver>, debug: bool) -> FeosResult<()> {
29        // Solve the profile
30        self.profile.solve(solver, debug)?;
31
32        // postprocess
33        self.surface_tension = Some(
34            (self.profile.integrate(
35                &(self.profile.grand_potential_density()?
36                    + self.vle.vapor().pressure(Contributions::Total)),
37            )) / Area::from_reduced(1.0),
38        );
39        let delta_rho = self.vle.liquid().density - self.vle.vapor().density;
40        self.equimolar_radius = Some(
41            self.profile
42                .integrate(&(self.profile.density.sum_axis(Axis_nd(0)) - self.vle.vapor().density))
43                / delta_rho
44                / Area::from_reduced(1.0),
45        );
46
47        Ok(())
48    }
49
50    pub fn solve(mut self, solver: Option<&DFTSolver>) -> FeosResult<Self> {
51        self.solve_inplace(solver, false)?;
52        Ok(self)
53    }
54}
55
56impl<F: HelmholtzEnergyFunctional> PlanarInterface<F> {
57    pub fn new(vle: &PhaseEquilibrium<F, 2>, n_grid: usize, l_grid: Length) -> Self {
58        // generate grid
59        let grid = Grid::Cartesian1(Axis::new_cartesian(n_grid, l_grid, None));
60
61        Self {
62            profile: DFTProfile::new(grid, vle.vapor(), None, None, None),
63            vle: vle.clone(),
64            surface_tension: None,
65            equimolar_radius: None,
66        }
67    }
68
69    pub fn from_tanh(
70        vle: &PhaseEquilibrium<F, 2>,
71        n_grid: usize,
72        l_grid: Length,
73        critical_temperature: Temperature,
74        fix_equimolar_surface: bool,
75    ) -> Self {
76        let mut profile = Self::new(vle, n_grid, l_grid);
77
78        // calculate segment indices
79        let indices = &profile.profile.bulk.eos.component_index();
80
81        // calculate density profile
82        let z0 = 0.5 * l_grid.to_reduced();
83        let (z0, sign) = (z0.abs(), -z0.signum());
84        let reduced_temperature = (vle.vapor().temperature / critical_temperature).into_value();
85        profile.profile.density =
86            Density::from_shape_fn(profile.profile.density.raw_dim(), |(i, z)| {
87                let rho_v = profile.vle.vapor().partial_density.get(indices[i]);
88                let rho_l = profile.vle.liquid().partial_density.get(indices[i]);
89                0.5 * (rho_l - rho_v)
90                    * (sign * (profile.profile.grid.grids()[0][z] - z0) / 3.0
91                        * (2.4728 - 2.3625 * reduced_temperature))
92                        .tanh()
93                    + 0.5 * (rho_l + rho_v)
94            });
95
96        // specify specification
97        if fix_equimolar_surface {
98            profile.profile.specification = Arc::new(DFTSpecifications::total_moles_from_profile(
99                &profile.profile,
100            ));
101        }
102
103        profile
104    }
105
106    pub fn from_pdgt(
107        vle: &PhaseEquilibrium<F, 2>,
108        n_grid: usize,
109        fix_equimolar_surface: bool,
110    ) -> FeosResult<Self> {
111        let dft = &vle.vapor().eos;
112
113        if dft.component_index().len() != 1 {
114            panic!("Initialization from pDGT not possible for segment DFT or mixtures");
115        }
116
117        // calculate density profile from pDGT
118        let n_grid_pdgt = 20;
119        let mut z_pdgt = Length::zeros(n_grid_pdgt);
120        let mut w_pdgt = Length::from_reduced(0.0);
121        let (rho_pdgt, gamma_pdgt) =
122            dft.solve_pdgt(vle, 20, 0, Some((&mut z_pdgt, &mut w_pdgt)))?;
123        if !gamma_pdgt.to_reduced().is_normal() {
124            return Err(FeosError::InvalidState(
125                String::from("DFTProfile::from_pdgt"),
126                String::from("gamma_pdgt"),
127                gamma_pdgt.to_reduced(),
128            ));
129        }
130
131        // create PlanarInterface
132        let l_grid = Length::from_reduced(MIN_WIDTH).max(w_pdgt * RELATIVE_WIDTH);
133        let mut profile = Self::new(vle, n_grid, l_grid);
134
135        // interpolate density profile from pDGT to DFT
136        let r = l_grid * 0.5;
137        profile.profile.density = interp_symmetric(
138            vle,
139            z_pdgt,
140            rho_pdgt,
141            &profile.vle,
142            profile.profile.grid.grids()[0],
143            r,
144        )?;
145
146        // specify specification
147        if fix_equimolar_surface {
148            profile.profile.specification = Arc::new(DFTSpecifications::total_moles_from_profile(
149                &profile.profile,
150            ));
151        }
152
153        Ok(profile)
154    }
155}
156
157impl<F: HelmholtzEnergyFunctional> PlanarInterface<F> {
158    pub fn shift_equimolar_inplace(&mut self) {
159        let s = self.profile.density.shape();
160        let m = &self.profile.bulk.eos.m();
161        let mut rho_l = Density::from_reduced(0.0);
162        let mut rho_v = Density::from_reduced(0.0);
163        let mut rho = Density::zeros(s[1]);
164        for i in 0..s[0] {
165            rho_l += self.profile.density.get((i, 0)) * m[i];
166            rho_v += self.profile.density.get((i, s[1] - 1)) * m[i];
167            rho += &(&self.profile.density.index_axis(Axis_nd(0), i) * m[i]);
168        }
169
170        let x = (rho - rho_v) / (rho_l - rho_v);
171        let ze = self.profile.grid.axes()[0].edges[0] + self.profile.integrate(&x).to_reduced();
172        self.profile.grid.axes_mut()[0].grid -= ze;
173    }
174
175    pub fn shift_equimolar(mut self) -> Self {
176        self.shift_equimolar_inplace();
177        self
178    }
179
180    /// Relative adsorption of component `i' with respect to `j': \Gamma_i^(j)
181    pub fn relative_adsorption(&self) -> Moles<Array2<f64>> {
182        let s = self.profile.density.shape();
183        let mut rho_l = Density::zeros(s[0]);
184        let mut rho_v = Density::zeros(s[0]);
185
186        // Calculate the partial densities in the liquid and in the vapor phase
187        for i in 0..s[0] {
188            rho_l.set(i, self.profile.density.get((i, 0)));
189            rho_v.set(i, self.profile.density.get((i, s[1] - 1)));
190        }
191
192        // Calculate \Gamma_i^(j)
193        Moles::from_shape_fn((s[0], s[0]), |(i, j)| {
194            if i == j {
195                Moles::from_reduced(0.0)
196            } else {
197                self.profile.integrate(
198                    &(-(rho_l.get(i) - rho_v.get(i))
199                        * ((&self.profile.density.index_axis(Axis_nd(0), j) - rho_l.get(j))
200                            / (rho_l.get(j) - rho_v.get(j))
201                            - (&self.profile.density.index_axis(Axis_nd(0), i) - rho_l.get(i))
202                                / (rho_l.get(i) - rho_v.get(i)))),
203                )
204            }
205        })
206    }
207
208    /// Interfacial enrichment of component `i': E_i
209    pub fn interfacial_enrichment(&self) -> Array1<f64> {
210        let s = self.profile.density.shape();
211        let density = self.profile.density.to_reduced();
212        let rho_l = density.index_axis(Axis_nd(1), 0);
213        let rho_v = density.index_axis(Axis_nd(1), s[1] - 1);
214
215        Array1::from_shape_fn(s[0], |i| {
216            *(density
217                .index_axis(Axis_nd(0), i)
218                .iter()
219                .max_by(|&a, &b| a.total_cmp(b))
220                .unwrap())  // panics only of iterator is empty
221                / rho_l[i].max(rho_v[i])
222        })
223    }
224
225    /// Interface thickness (90-10 number density difference)
226    pub fn interfacial_thickness(&self) -> FeosResult<Length> {
227        let s = self.profile.density.shape();
228        let rho = self.profile.density.sum_axis(Axis_nd(0)).to_reduced();
229        let z = self.profile.grid.grids()[0];
230        let dz = z[1] - z[0];
231
232        let limits = (0.9_f64, 0.1_f64);
233        let (limit_upper, limit_lower) = if limits.0 > limits.1 {
234            (limits.0, limits.1)
235        } else {
236            (limits.1, limits.0)
237        };
238
239        if limit_upper >= 1.0 || limit_upper.is_sign_negative() {
240            return Err(FeosError::IterationFailed(String::from(
241                "Upper limit 'l' of interface thickness needs to satisfy 0 < l < 1.",
242            )));
243        }
244        if limit_lower >= 1.0 || limit_lower.is_sign_negative() {
245            return Err(FeosError::IterationFailed(String::from(
246                "Lower limit 'l' of interface thickness needs to satisfy 0 < l < 1.",
247            )));
248        }
249
250        // Get the densities in the liquid and in the vapor phase
251        let rho_v = rho[0].min(rho[s[1] - 1]);
252        let rho_l = rho[0].max(rho[s[1] - 1]);
253
254        if (rho_l - rho_v).abs() < 1.0e-10 {
255            return Ok(Length::from_reduced(0.0));
256        }
257
258        // Density boundaries for interface definition
259        let rho_upper = rho_v + limit_upper * (rho_l - rho_v);
260        let rho_lower = rho_v + limit_lower * (rho_l - rho_v);
261
262        // Get indizes right of intersection between density profile and
263        // constant density boundaries
264        let index_upper_plus = if rho[0] >= rho[s[1] - 1] {
265            rho.iter()
266                .enumerate()
267                .find(|&(_, &x)| (x - rho_upper).is_sign_negative())
268                .expect("Could not find rho_upper value!")
269                .0
270        } else {
271            rho.iter()
272                .enumerate()
273                .find(|&(_, &x)| (rho_upper - x).is_sign_negative())
274                .expect("Could not find rho_upper value!")
275                .0
276        };
277        let index_lower_plus = if rho[0] >= rho[s[1] - 1] {
278            rho.iter()
279                .enumerate()
280                .find(|&(_, &x)| (x - rho_lower).is_sign_negative())
281                .expect("Could not find rho_lower value!")
282                .0
283        } else {
284            rho.iter()
285                .enumerate()
286                .find(|&(_, &x)| (rho_lower - x).is_sign_negative())
287                .expect("Could not find rho_lower value!")
288                .0
289        };
290
291        // Calculate distance between two density points using a linear
292        // interpolated density profiles between the two grid points where the
293        // density profile crosses the limiting densities
294        let z_upper = z[index_upper_plus - 1]
295            + (rho_upper - rho[index_upper_plus - 1])
296                / (rho[index_upper_plus] - rho[index_upper_plus - 1])
297                * dz;
298        let z_lower = z[index_lower_plus - 1]
299            + (rho_lower - rho[index_lower_plus - 1])
300                / (rho[index_lower_plus] - rho[index_lower_plus - 1])
301                * dz;
302
303        // Return
304        Ok(Length::from_reduced(z_lower - z_upper))
305    }
306
307    fn set_density_scale(&mut self, init: &Density<Array2<f64>>) {
308        assert_eq!(self.profile.density.shape(), init.shape());
309        let n_grid = self.profile.density.shape()[1];
310        let drho_init = &init.index_axis(Axis_nd(1), 0) - &init.index_axis(Axis_nd(1), n_grid - 1);
311        let rho_init_0 = init.index_axis(Axis_nd(1), n_grid - 1);
312        let drho = &self.profile.density.index_axis(Axis_nd(1), 0)
313            - &self.profile.density.index_axis(Axis_nd(1), n_grid - 1);
314        let rho_0 = self.profile.density.index_axis(Axis_nd(1), n_grid - 1);
315
316        self.profile.density = Density::from_shape_fn(self.profile.density.raw_dim(), |(i, j)| {
317            ((init.get((i, j)) - rho_init_0.get(i)) / drho_init.get(i)).into_value() * drho.get(i)
318                + rho_0.get(i)
319        });
320    }
321
322    pub fn set_density_inplace(&mut self, init: &Density<Array2<f64>>, scale: bool) {
323        if scale {
324            self.set_density_scale(init)
325        } else {
326            assert_eq!(self.profile.density.shape(), init.shape());
327            self.profile.density = init.clone();
328        }
329    }
330
331    pub fn set_density(mut self, init: &Density<Array2<f64>>, scale: bool) -> Self {
332        self.set_density_inplace(init, scale);
333        self
334    }
335}
336
337fn interp_symmetric<F: HelmholtzEnergyFunctional>(
338    vle_pdgt: &PhaseEquilibrium<F, 2>,
339    z_pdgt: Length<Array1<f64>>,
340    rho_pdgt: Density<Array2<f64>>,
341    vle: &PhaseEquilibrium<F, 2>,
342    z: &Array1<f64>,
343    radius: Length,
344) -> FeosResult<Density<Array2<f64>>> {
345    let reduced_density = Array2::from_shape_fn(rho_pdgt.raw_dim(), |(i, j)| {
346        ((rho_pdgt.get((i, j)) - vle_pdgt.vapor().partial_density.get(i))
347            / (vle_pdgt.liquid().partial_density.get(i) - vle_pdgt.vapor().partial_density.get(i)))
348        .into_value()
349            - 0.5
350    });
351    let segments = vle_pdgt.vapor().eos.component_index().len();
352    let mut reduced_density = interp(
353        &z_pdgt.to_reduced(),
354        &reduced_density,
355        &(z - radius.to_reduced()),
356        &Array1::from_elem(segments, 0.5),
357        &Array1::from_elem(segments, -0.5),
358        false,
359    ) + interp(
360        &z_pdgt.to_reduced(),
361        &reduced_density,
362        &(z + radius.to_reduced()),
363        &Array1::from_elem(segments, -0.5),
364        &Array1::from_elem(segments, 0.5),
365        true,
366    );
367    if radius.is_sign_negative() {
368        reduced_density += 1.0;
369    }
370    Ok(Density::from_shape_fn(
371        reduced_density.raw_dim(),
372        |(i, j)| {
373            reduced_density[(i, j)]
374                * (vle.liquid().partial_density.get(i) - vle.vapor().partial_density.get(i))
375                + vle.vapor().partial_density.get(i)
376        },
377    ))
378}
379
380fn interp(
381    x_old: &Array1<f64>,
382    y_old: &Array2<f64>,
383    x_new: &Array1<f64>,
384    y_left: &Array1<f64>,
385    y_right: &Array1<f64>,
386    reverse: bool,
387) -> Array2<f64> {
388    let n = x_old.len();
389
390    let (x_rev, y_rev) = if reverse {
391        (-&x_old.slice(s![..;-1]), y_old.slice(s![.., ..;-1]))
392    } else {
393        (x_old.to_owned(), y_old.view())
394    };
395
396    let mut y_new = Array2::zeros((y_rev.shape()[0], x_new.len()));
397    let mut k = 0;
398    for i in 0..x_new.len() {
399        while k < n && x_new[i] > x_rev[k] {
400            k += 1;
401        }
402        y_new.slice_mut(s![.., i]).assign(&if k == 0 {
403            y_left
404                + &((&y_rev.slice(s![.., 0]) - y_left)
405                    * ((&y_rev.slice(s![.., 1]) - y_left) / (&y_rev.slice(s![.., 0]) - y_left))
406                        .mapv(|x| x.powf((x_new[i] - x_rev[0]) / (x_rev[1] - x_rev[0]))))
407        } else if k == n {
408            y_right
409                + &((&y_rev.slice(s![.., n - 2]) - y_right)
410                    * ((&y_rev.slice(s![.., n - 1]) - y_right)
411                        / (&y_rev.slice(s![.., n - 2]) - y_right))
412                        .mapv(|x| {
413                            x.powf((x_new[i] - x_rev[n - 2]) / (x_rev[n - 1] - x_rev[n - 2]))
414                        }))
415        } else {
416            &y_rev.slice(s![.., k - 1])
417                + &((x_new[i] - x_rev[k - 1]) / (x_rev[k] - x_rev[k - 1])
418                    * (&y_rev.slice(s![.., k]) - &y_rev.slice(s![.., k - 1])))
419        });
420    }
421    y_new
422}
feos_dft/interface/mod.rs

feos_dft/interface/
mod.rs
