Skip to content

Commit 447dae1

Browse files
prehnerprehner
committed
cleanup and changelog entry
1 parent af247ac commit 447dae1

File tree

2 files changed

+25
-26
lines changed

2 files changed

+25
-26
lines changed

‎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
## [Breaking]
8+
### Added
9+
- Extended tp-flash algorithm to static numbers of components and enabled automatic differentiation for binary systems. [#336](https://github.com/feos-org/feos/pull/336)
10+
811
### Packaging
912
- Updated `quantity` dependency to 0.13 and removed the `typenum` dependency. [#323](https://github.com/feos-org/feos/pull/323)
1013

‎crates/feos-core/src/phase_equilibria/tp_flash.rs‎

Lines changed: 22 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,9 @@ use crate::state::{Contributions, State};
55
use crate::{ReferenceSystem, SolverOptions, Verbosity};
66
use nalgebra::allocator::Allocator;
77
use nalgebra::{DefaultAllocator, Dim, Matrix3, OVector, SVector, U1, U2, vector};
8-
use num_dual::{Dual, DualNum, DualStruct, Gradients, first_derivative, implicit_derivative_sp};
8+
use num_dual::{
9+
Dual, Dual2Vec, DualNum, DualStruct, Gradients, first_derivative, implicit_derivative_sp,
10+
};
911
use quantity::{Dimensionless, MOL, MolarVolume, Moles, Pressure, Quantity, Temperature};
1012

1113
const MAX_ITER_TP: usize = 400;
@@ -39,11 +41,10 @@ where
3941
}
4042

4143
impl<E: Residual<U2, D>, D: DualNum<f64> + Copy> PhaseEquilibrium<E, 2, U2, D> {
42-
/// Perform a Tp-flash calculation. If no initial values are
43-
/// given, the solution is initialized using a stability analysis.
44-
///
45-
/// The algorithm can be use to calculate phase equilibria of systems
46-
/// containing non-volatile components (e.g. ions).
44+
/// Perform a Tp-flash calculation for a binary mixture.
45+
/// Compared to the version of the algorithm for a generic
46+
/// number of components ([tp_flash](PhaseEquilibrium::tp_flash)),
47+
/// this can be used in combination with automatic differentiation.
4748
pub fn tp_flash_binary(
4849
eos: &E,
4950
temperature: Temperature<D>,
@@ -65,34 +66,29 @@ impl<E: Residual<U2, D>, D: DualNum<f64> + Copy> PhaseEquilibrium<E, 2, U2, D> {
6566
let x = vle_re.liquid().molefracs[0];
6667
let y = vle_re.vapor().molefracs[0];
6768
let x = implicit_derivative_sp(
68-
|x, args: &SVector<_, _>| {
69+
|x, &[t, p, z]: &[_; 3]| {
6970
let [[v_l, v_v, x, y]] = x.data.0;
70-
let [[t, p, z]] = args.data.0;
7171
let beta = (z - x) / (y - x);
72-
let molefracs_l = vector![x, -x + 1.0];
73-
let molefracs_v = vector![y, -y + 1.0];
7472
let eos = eos.lift();
75-
let a_l_res = eos.residual_molar_helmholtz_energy(t, v_l, &molefracs_l);
76-
let a_l_ig = (x * (x / v_l).ln() - (x - 1.0) * ((-x + 1.0) / v_l).ln() - 1.0) * t;
77-
let a_v_res = eos.residual_molar_helmholtz_energy(t, v_v, &molefracs_v);
78-
let a_v_ig = (y * (y / v_v).ln() - (y - 1.0) * ((-y + 1.0) / v_v).ln() - 1.0) * t;
79-
let a_res = (a_v_res + a_v_ig) * beta - (a_l_res + a_l_ig) * (beta - 1.0);
80-
let v = v_v * beta - v_l * (beta - 1.0);
81-
a_res + v * p
73+
let pot = |x: Dual2Vec<_, _, _>, v| {
74+
let molefracs = vector![x, -x + 1.0];
75+
let a_res = eos.residual_molar_helmholtz_energy(t, v, &molefracs);
76+
let a_ig = (x * (x / v).ln() - (x - 1.0) * ((-x + 1.0) / v).ln() - 1.0) * t;
77+
a_res + a_ig + v * p
78+
};
79+
pot(y, v_v) * beta - pot(x, v_l) * (beta - 1.0)
8280
},
8381
SVector::from([v_l, v_v, x, y]),
84-
&SVector::from([t, p, z]),
82+
&[t, p, z],
8583
);
8684
let [[v_l, v_v, x, y]] = x.data.0;
8785
let beta = (z - x) / (y - x);
88-
let volume = MolarVolume::from_reduced(v_l * (-beta + 1.0)) * total_moles;
89-
let moles =
90-
Quantity::new(vector![x, -x + 1.0] * (-beta + 1.0) * total_moles.convert_into(MOL));
91-
let liquid = State::new_nvt(eos, temperature, volume, &moles)?;
92-
let volume = MolarVolume::from_reduced(v_v * beta) * total_moles;
93-
let moles = Quantity::new(vector![y, -y + 1.0] * beta * total_moles.convert_into(MOL));
94-
let vapor = State::new_nvt(eos, temperature, volume, &moles)?;
95-
Ok(Self([vapor, liquid]))
86+
let state = |x: D, v, phi| {
87+
let volume = MolarVolume::from_reduced(v * phi) * total_moles;
88+
let moles = Quantity::new(vector![x, -x + 1.0] * phi * total_moles.convert_into(MOL));
89+
State::new_nvt(eos, temperature, volume, &moles)
90+
};
91+
Ok(Self([state(y, v_v, beta)?, state(x, v_l, -beta + 1.0)?]))
9692
}
9793
}
9894

0 commit comments

Comments
 (0)