Skip to content

Commit fc19a14

Browse files
prehnerprehner
committed
Implement ph and ps flashes for binary mixtures
1 parent c33df59 commit fc19a14

File tree

16 files changed

+1082
-94
lines changed

16 files changed

+1082
-94
lines changed

‎CHANGELOG.md‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -10,15 +10,18 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
1010
- Rewrote `PhaseEquilibrium::pure_p` to mirror `pure_t` and enabled automatic differentiation. [#337](https://github.com/feos-org/feos/pull/337)
1111
- Added `boiling_temperature` to the list of properties for parallel evaluations of gradients. [#337](https://github.com/feos-org/feos/pull/337)
1212
- Added the `Composition` trait to allow more flexibility in the creation of states and phase equilibria. [#330](https://github.com/feos-org/feos/pull/330)
13+
- Added `PhaseEquilibrium::ph_flash` and `PhaseEquilibrium::ps_flash`. [#338](https://github.com/feos-org/feos/pull/338)
14+
- Added getters for `vapor_phase_fraction`, `molar_enthalpy`, `molar_entropy`, `total_moles`, `enthalpy`, and `entropy` to `PhaseEquilibrium`. [#338](https://github.com/feos-org/feos/pull/338)
1315

1416
### Changed
1517
- Removed any assumptions about the total number of moles in a `State` or `PhaseEquilibrium`. Evaluating extensive properties now returns a `Result`. [#330](https://github.com/feos-org/feos/pull/330)
18+
- Redesigned the `IdealGas` trait and added `IdealGasAD` in analogy to `ResidualDyn` and `Residual`. [#330](https://github.com/feos-org/feos/pull/330)
1619

1720
### Removed
1821
- Removed the `StateBuilder` struct, because it is mostly obsolete with the addition of the `Composition` trait. [#330](https://github.com/feos-org/feos/pull/330)
1922

2023
### Packaging
21-
- Updated `quantity` dependency to 0.13 and removed the `typenum` dependency. [#323](https://github.com/feos-org/feos/pull/323)
24+
- Updated `quantity` dependency to 0.13 and removed the `typenum` dependency. [#328](https://github.com/feos-org/feos/pull/328)
2225

2326
## [Unreleased]
2427
### Added

‎crates/feos-core/src/equation_of_state/mod.rs‎

Lines changed: 79 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
use crate::ReferenceSystem;
22
use crate::state::StateHD;
33
use nalgebra::{
4-
Const, DVector, DefaultAllocator, Dim, Dyn, OVector, SVector, U1, allocator::Allocator,
4+
Const, DVector, DefaultAllocator, Dim, Dyn, OVector, SVector, allocator::Allocator,
55
};
66
use num_dual::DualNum;
77
use quantity::{Dimensionless, MolarEnergy, MolarVolume, Temperature};
@@ -114,11 +114,28 @@ impl<I: Clone, R: Residual<Const<N>, D>, D: DualNum<f64> + Copy, const N: usize>
114114
}
115115

116116
/// Ideal gas Helmholtz energy contribution.
117-
pub trait IdealGas<D = f64> {
117+
pub trait IdealGas {
118118
/// Implementation of an ideal gas model in terms of the
119119
/// logarithm of the cubic thermal de Broglie wavelength
120120
/// in units ln(A³) for each component in the system.
121-
fn ln_lambda3<D2: DualNum<f64, Inner = D> + Copy>(&self, temperature: D2) -> D2;
121+
fn ln_lambda3<D: DualNum<f64> + Copy>(&self, temperature: D) -> D;
122+
123+
/// The name of the ideal gas model.
124+
fn ideal_gas_model(&self) -> &'static str;
125+
}
126+
127+
/// Ideal gas Helmholtz energy contribution with automatic differentiation with
128+
/// respect to parameters.
129+
pub trait IdealGasAD<D = f64>: Clone {
130+
type Real: IdealGasAD;
131+
type Lifted<D2: DualNum<f64, Inner = D> + Copy>: IdealGasAD<D2>;
132+
fn re(&self) -> Self::Real;
133+
fn lift<D2: DualNum<f64, Inner = D> + Copy>(&self) -> Self::Lifted<D2>;
134+
135+
/// Implementation of an ideal gas model in terms of the
136+
/// logarithm of the cubic thermal de Broglie wavelength
137+
/// in units ln(A³) for each component in the system.
138+
fn ln_lambda3(&self, temperature: D) -> D;
122139

123140
/// The name of the ideal gas model.
124141
fn ideal_gas_model(&self) -> &'static str;
@@ -129,45 +146,44 @@ pub trait Total<N: Dim = Dyn, D: DualNum<f64> + Copy = f64>: Residual<N, D>
129146
where
130147
DefaultAllocator: Allocator<N>,
131148
{
132-
type IdealGas: IdealGas<D>;
149+
type RealTotal: Total<N, f64>;
150+
type LiftedTotal<D2: DualNum<f64, Inner = D> + Copy>: Total<N, D2>;
151+
fn re_total(&self) -> Self::RealTotal;
152+
fn lift_total<D2: DualNum<f64, Inner = D> + Copy>(&self) -> Self::LiftedTotal<D2>;
133153

134154
fn ideal_gas_model(&self) -> &'static str;
135155

136-
fn ideal_gas(&self) -> impl Iterator<Item = &Self::IdealGas>;
137-
138-
fn ln_lambda3<D2: DualNum<f64, Inner = D> + Copy>(&self, temperature: D2) -> OVector<D2, N> {
139-
OVector::from_iterator_generic(
140-
N::from_usize(self.components()),
141-
U1,
142-
self.ideal_gas().map(|i| i.ln_lambda3(temperature)),
143-
)
144-
}
156+
fn ln_lambda3(&self, temperature: D) -> OVector<D, N>;
145157

146-
fn ideal_gas_molar_helmholtz_energy<D2: DualNum<f64, Inner = D> + Copy>(
158+
fn ideal_gas_molar_helmholtz_energy(
147159
&self,
148-
temperature: D2,
149-
molar_volume: D2,
150-
molefracs: &OVector<D2, N>,
151-
) -> D2 {
160+
temperature: D,
161+
molar_volume: D,
162+
molefracs: &OVector<D, N>,
163+
) -> D {
152164
let partial_density = molefracs / molar_volume;
153-
let mut res = D2::from(0.0);
154-
for (i, &r) in self.ideal_gas().zip(partial_density.iter()) {
165+
let mut res = D::from(0.0);
166+
for (&l, &r) in self
167+
.ln_lambda3(temperature)
168+
.iter()
169+
.zip(partial_density.iter())
170+
{
155171
let ln_rho_m1 = if r.re() == 0.0 {
156-
D2::from(0.0)
172+
D::from(0.0)
157173
} else {
158174
r.ln() - 1.0
159175
};
160-
res += r * (i.ln_lambda3(temperature) + ln_rho_m1)
176+
res += r * (l + ln_rho_m1)
161177
}
162178
res * molar_volume * temperature
163179
}
164180

165-
fn ideal_gas_helmholtz_energy<D2: DualNum<f64, Inner = D> + Copy>(
181+
fn ideal_gas_helmholtz_energy(
166182
&self,
167-
temperature: Temperature<D2>,
168-
volume: MolarVolume<D2>,
169-
moles: &OVector<D2, N>,
170-
) -> MolarEnergy<D2> {
183+
temperature: Temperature<D>,
184+
volume: MolarVolume<D>,
185+
moles: &OVector<D, N>,
186+
) -> MolarEnergy<D> {
171187
let total_moles = moles.sum();
172188
let molefracs = moles / total_moles;
173189
let molar_volume = volume.into_reduced() / total_moles;
@@ -180,32 +196,59 @@ where
180196
}
181197

182198
impl<
183-
I: IdealGas + Clone + 'static,
199+
I: IdealGas + 'static,
184200
C: Deref<Target = EquationOfState<Vec<I>, R>> + Clone,
185201
R: ResidualDyn + 'static,
186-
> Total<Dyn, f64> for C
202+
D: DualNum<f64> + Copy,
203+
> Total<Dyn, D> for C
187204
{
188-
type IdealGas = I;
205+
type RealTotal = Self;
206+
type LiftedTotal<D2: DualNum<f64, Inner = D> + Copy> = Self;
207+
fn re_total(&self) -> Self::RealTotal {
208+
self.clone()
209+
}
210+
fn lift_total<D2: DualNum<f64, Inner = D> + Copy>(&self) -> Self::LiftedTotal<D2> {
211+
self.clone()
212+
}
189213

190214
fn ideal_gas_model(&self) -> &'static str {
191215
self.ideal_gas[0].ideal_gas_model()
192216
}
193217

194-
fn ideal_gas(&self) -> impl Iterator<Item = &Self::IdealGas> {
195-
self.ideal_gas.iter()
218+
fn ln_lambda3(&self, temperature: D) -> DVector<D> {
219+
DVector::from_vec(
220+
self.ideal_gas
221+
.iter()
222+
.map(|i| i.ln_lambda3(temperature))
223+
.collect(),
224+
)
196225
}
197226
}
198227

199-
impl<I: IdealGas<D> + Clone, R: Residual<Const<N>, D>, D: DualNum<f64> + Copy, const N: usize>
228+
impl<I: IdealGasAD<D>, R: Residual<Const<N>, D>, D: DualNum<f64> + Copy, const N: usize>
200229
Total<Const<N>, D> for EquationOfState<[I; N], R>
201230
{
202-
type IdealGas = I;
231+
type RealTotal = EquationOfState<[I::Real; N], R::Real>;
232+
type LiftedTotal<D2: DualNum<f64, Inner = D> + Copy> =
233+
EquationOfState<[I::Lifted<D2>; N], R::Lifted<D2>>;
234+
fn re_total(&self) -> Self::RealTotal {
235+
EquationOfState::new(
236+
self.ideal_gas.each_ref().map(|i| i.re()),
237+
self.residual.re(),
238+
)
239+
}
240+
fn lift_total<D2: DualNum<f64, Inner = D> + Copy>(&self) -> Self::LiftedTotal<D2> {
241+
EquationOfState::new(
242+
self.ideal_gas.each_ref().map(|i| i.lift()),
243+
self.residual.lift(),
244+
)
245+
}
203246

204247
fn ideal_gas_model(&self) -> &'static str {
205248
self.ideal_gas[0].ideal_gas_model()
206249
}
207250

208-
fn ideal_gas(&self) -> impl Iterator<Item = &Self::IdealGas> {
209-
self.ideal_gas.iter()
251+
fn ln_lambda3(&self, temperature: D) -> SVector<D, N> {
252+
SVector::from(self.ideal_gas.each_ref().map(|i| i.ln_lambda3(temperature)))
210253
}
211254
}

‎crates/feos-core/src/lib.rs‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -34,8 +34,8 @@ mod phase_equilibria;
3434
mod state;
3535
pub use ad::{ParametersAD, PropertiesAD};
3636
pub use equation_of_state::{
37-
EntropyScaling, EquationOfState, IdealGas, Molarweight, NoResidual, Residual, ResidualDyn,
38-
Subset, Total,
37+
EntropyScaling, EquationOfState, IdealGas, IdealGasAD, Molarweight, NoResidual, Residual,
38+
ResidualDyn, Subset, Total,
3939
};
4040
pub use errors::{FeosError, FeosResult};
4141
#[cfg(feature = "ndarray")]

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

Lines changed: 27 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -10,16 +10,23 @@ use quantity::{Dimensionless, Energy, Entropy, MolarEnergy, MolarEntropy, Moles}
1010
use std::fmt;
1111
use std::fmt::Write;
1212

13+
// with empty lines to not mess up the order in the documentation
14+
mod vle_pure;
15+
1316
mod bubble_dew;
17+
18+
mod tp_flash;
19+
20+
mod px_flashes;
21+
1422
#[cfg(feature = "ndarray")]
1523
mod phase_diagram_binary;
1624
#[cfg(feature = "ndarray")]
1725
mod phase_diagram_pure;
1826
#[cfg(feature = "ndarray")]
1927
mod phase_envelope;
2028
mod stability_analysis;
21-
mod tp_flash;
22-
mod vle_pure;
29+
2330
pub use bubble_dew::TemperatureOrPressure;
2431
#[cfg(feature = "ndarray")]
2532
pub use phase_diagram_binary::PhaseDiagramHetero;
@@ -33,22 +40,25 @@ pub use phase_diagram_pure::PhaseDiagram;
3340
///
3441
/// ## Contents
3542
///
43+
/// + [Pure component phase equilibria](#pure-component-phase-equilibria)
3644
/// + [Bubble and dew point calculations](#bubble-and-dew-point-calculations)
37-
/// + [Heteroazeotropes](#heteroazeotropes)
3845
/// + [Flash calculations](#flash-calculations)
39-
/// + [Pure component phase equilibria](#pure-component-phase-equilibria)
46+
/// + [Heteroazeotropes](#heteroazeotropes)
4047
/// + [Utility functions](#utility-functions)
4148
#[derive(Debug, Clone)]
4249
pub struct PhaseEquilibrium<E, const P: usize, N: Dim = Dyn, D: DualNum<f64> + Copy = f64>
4350
where
4451
DefaultAllocator: Allocator<N>,
4552
{
46-
states: [State<E, N, D>; P],
53+
pub states: [State<E, N, D>; P],
4754
pub phase_fractions: [D; P],
4855
total_moles: Option<Moles<D>>,
4956
}
5057

51-
impl<E: Residual, const P: usize> fmt::Display for PhaseEquilibrium<E, P> {
58+
impl<E: Residual<N>, N: Dim, const P: usize> fmt::Display for PhaseEquilibrium<E, P, N>
59+
where
60+
DefaultAllocator: Allocator<N>,
61+
{
5262
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
5363
for (i, s) in self.states.iter().enumerate() {
5464
writeln!(f, "phase {i}: {s}")?;
@@ -125,12 +135,12 @@ impl<E: Residual<N, D>, N: Dim, D: DualNum<f64> + Copy> PhaseEquilibrium<E, 2, N
125135
where
126136
DefaultAllocator: Allocator<N>,
127137
{
128-
pub(super) fn single_phase(state: State<E, N, D>) -> Self {
138+
pub fn single_phase(state: State<E, N, D>) -> Self {
129139
let total_moles = state.total_moles;
130140
Self::with_vapor_phase_fraction(state.clone(), state, D::from(1.0), total_moles)
131141
}
132142

133-
pub(super) fn two_phase(vapor: State<E, N, D>, liquid: State<E, N, D>) -> Self {
143+
pub fn two_phase(vapor: State<E, N, D>, liquid: State<E, N, D>) -> Self {
134144
let (beta, total_moles) =
135145
if let (Some(nv), Some(nl)) = (vapor.total_moles, liquid.total_moles) {
136146
(nv.convert_into(nl + nv), Some(nl + nv))
@@ -140,7 +150,7 @@ where
140150
Self::with_vapor_phase_fraction(vapor, liquid, beta, total_moles)
141151
}
142152

143-
pub(super) fn with_vapor_phase_fraction(
153+
pub fn with_vapor_phase_fraction(
144154
vapor: State<E, N, D>,
145155
liquid: State<E, N, D>,
146156
vapor_phase_fraction: D,
@@ -158,11 +168,7 @@ impl<E: Residual<N, D>, N: Dim, D: DualNum<f64> + Copy> PhaseEquilibrium<E, 3, N
158168
where
159169
DefaultAllocator: Allocator<N>,
160170
{
161-
pub(super) fn new(
162-
vapor: State<E, N, D>,
163-
liquid1: State<E, N, D>,
164-
liquid2: State<E, N, D>,
165-
) -> Self {
171+
pub fn new(vapor: State<E, N, D>, liquid1: State<E, N, D>, liquid2: State<E, N, D>) -> Self {
166172
Self {
167173
states: [vapor, liquid1, liquid2],
168174
phase_fractions: [D::from(1.0), D::from(0.0), D::from(0.0)],
@@ -171,15 +177,21 @@ where
171177
}
172178
}
173179

174-
impl<E: Total<N, D>, N: Gradients, const P: usize, D: DualNum<f64> + Copy>
180+
impl<E: Residual<N, D>, N: Gradients, const P: usize, D: DualNum<f64> + Copy>
175181
PhaseEquilibrium<E, P, N, D>
176182
where
177183
DefaultAllocator: Allocator<N>,
178184
{
179185
pub fn total_moles(&self) -> FeosResult<Moles<D>> {
180186
self.total_moles.ok_or(FeosError::IntensiveState)
181187
}
188+
}
182189

190+
impl<E: Total<N, D>, N: Gradients, const P: usize, D: DualNum<f64> + Copy>
191+
PhaseEquilibrium<E, P, N, D>
192+
where
193+
DefaultAllocator: Allocator<N>,
194+
{
183195
pub fn molar_enthalpy(&self) -> MolarEnergy<D> {
184196
self.states
185197
.iter()

0 commit comments

Comments
 (0)