Skip to content

Commit 8b52d56

Browse files
prehnerprehner
committed
Implement ph and ps flashes for binary mixtures
1 parent f609159 commit 8b52d56

File tree

12 files changed

+661
-69
lines changed

12 files changed

+661
-69
lines changed

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

Lines changed: 78 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,27 @@ 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.
128+
pub trait IdealGasAD<D = f64>: Clone {
129+
type Real: IdealGasAD;
130+
type Lifted<D2: DualNum<f64, Inner = D> + Copy>: IdealGasAD<D2>;
131+
fn re(&self) -> Self::Real;
132+
fn lift<D2: DualNum<f64, Inner = D> + Copy>(&self) -> Self::Lifted<D2>;
133+
134+
/// Implementation of an ideal gas model in terms of the
135+
/// logarithm of the cubic thermal de Broglie wavelength
136+
/// in units ln(A³) for each component in the system.
137+
fn ln_lambda3(&self, temperature: D) -> D;
122138

123139
/// The name of the ideal gas model.
124140
fn ideal_gas_model(&self) -> &'static str;
@@ -129,45 +145,44 @@ pub trait Total<N: Dim = Dyn, D: DualNum<f64> + Copy = f64>: Residual<N, D>
129145
where
130146
DefaultAllocator: Allocator<N>,
131147
{
132-
type IdealGas: IdealGas<D>;
148+
type RealTotal: Total<N, f64>;
149+
type LiftedTotal<D2: DualNum<f64, Inner = D> + Copy>: Total<N, D2>;
150+
fn re_total(&self) -> Self::RealTotal;
151+
fn lift_total<D2: DualNum<f64, Inner = D> + Copy>(&self) -> Self::LiftedTotal<D2>;
133152

134153
fn ideal_gas_model(&self) -> &'static str;
135154

136-
fn ideal_gas(&self) -> impl Iterator<Item = &Self::IdealGas>;
155+
fn ln_lambda3(&self, temperature: D) -> OVector<D, N>;
137156

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-
}
145-
146-
fn ideal_gas_molar_helmholtz_energy<D2: DualNum<f64, Inner = D> + Copy>(
157+
fn ideal_gas_molar_helmholtz_energy(
147158
&self,
148-
temperature: D2,
149-
molar_volume: D2,
150-
molefracs: &OVector<D2, N>,
151-
) -> D2 {
159+
temperature: D,
160+
molar_volume: D,
161+
molefracs: &OVector<D, N>,
162+
) -> D {
152163
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()) {
164+
let mut res = D::from(0.0);
165+
for (&l, &r) in self
166+
.ln_lambda3(temperature)
167+
.iter()
168+
.zip(partial_density.iter())
169+
{
155170
let ln_rho_m1 = if r.re() == 0.0 {
156-
D2::from(0.0)
171+
D::from(0.0)
157172
} else {
158173
r.ln() - 1.0
159174
};
160-
res += r * (i.ln_lambda3(temperature) + ln_rho_m1)
175+
res += r * (l + ln_rho_m1)
161176
}
162177
res * molar_volume * temperature
163178
}
164179

165-
fn ideal_gas_helmholtz_energy<D2: DualNum<f64, Inner = D> + Copy>(
180+
fn ideal_gas_helmholtz_energy(
166181
&self,
167-
temperature: Temperature<D2>,
168-
volume: MolarVolume<D2>,
169-
moles: &OVector<D2, N>,
170-
) -> MolarEnergy<D2> {
182+
temperature: Temperature<D>,
183+
volume: MolarVolume<D>,
184+
moles: &OVector<D, N>,
185+
) -> MolarEnergy<D> {
171186
let total_moles = moles.sum();
172187
let molefracs = moles / total_moles;
173188
let molar_volume = volume.into_reduced() / total_moles;
@@ -180,32 +195,59 @@ where
180195
}
181196

182197
impl<
183-
I: IdealGas + Clone + 'static,
198+
I: IdealGas + 'static,
184199
C: Deref<Target = EquationOfState<Vec<I>, R>> + Clone,
185200
R: ResidualDyn + 'static,
186-
> Total<Dyn, f64> for C
201+
D: DualNum<f64> + Copy,
202+
> Total<Dyn, D> for C
187203
{
188-
type IdealGas = I;
204+
type RealTotal = Self;
205+
type LiftedTotal<D2: DualNum<f64, Inner = D> + Copy> = Self;
206+
fn re_total(&self) -> Self::RealTotal {
207+
self.clone()
208+
}
209+
fn lift_total<D2: DualNum<f64, Inner = D> + Copy>(&self) -> Self::LiftedTotal<D2> {
210+
self.clone()
211+
}
189212

190213
fn ideal_gas_model(&self) -> &'static str {
191214
self.ideal_gas[0].ideal_gas_model()
192215
}
193216

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

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

204246
fn ideal_gas_model(&self) -> &'static str {
205247
self.ideal_gas[0].ideal_gas_model()
206248
}
207249

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

‎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: 22 additions & 6 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,10 +40,10 @@ 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>
@@ -48,7 +55,10 @@ where
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}")?;
@@ -171,15 +181,21 @@ where
171181
}
172182
}
173183

174-
impl<E: Total<N, D>, N: Gradients, const P: usize, D: DualNum<f64> + Copy>
184+
impl<E: Residual<N, D>, N: Gradients, const P: usize, D: DualNum<f64> + Copy>
175185
PhaseEquilibrium<E, P, N, D>
176186
where
177187
DefaultAllocator: Allocator<N>,
178188
{
179189
pub fn total_moles(&self) -> FeosResult<Moles<D>> {
180190
self.total_moles.ok_or(FeosError::IntensiveState)
181191
}
192+
}
182193

194+
impl<E: Total<N, D>, N: Gradients, const P: usize, D: DualNum<f64> + Copy>
195+
PhaseEquilibrium<E, P, N, D>
196+
where
197+
DefaultAllocator: Allocator<N>,
198+
{
183199
pub fn molar_enthalpy(&self) -> MolarEnergy<D> {
184200
self.states
185201
.iter()

0 commit comments

Comments
 (0)