Skip to content

Commit 6577305

Browse files
prehnerprehner
committed
include phase fraction in PhaseEquilibrium
1 parent 33ececb commit 6577305

32 files changed

+554
-608
lines changed

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

Lines changed: 13 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
use crate::DensityInitialization::Liquid;
22
use crate::density_iteration::density_iteration;
3-
use crate::{FeosResult, PhaseEquilibrium, ReferenceSystem, Residual};
3+
use crate::{Composition, FeosResult, PhaseEquilibrium, ReferenceSystem, Residual};
44
use nalgebra::{Const, SVector, U1, U2};
55
#[cfg(feature = "rayon")]
66
use ndarray::{Array1, Array2, ArrayView2, Zip};
@@ -215,20 +215,21 @@ pub trait PropertiesAD {
215215
)
216216
}
217217

218-
fn bubble_point_pressure<const P: usize>(
218+
fn bubble_point_pressure<const P: usize, X: Composition<f64, U2>>(
219219
&self,
220220
temperature: Temperature,
221221
pressure: Option<Pressure>,
222-
liquid_molefracs: SVector<f64, 2>,
222+
liquid_molefracs: X,
223223
) -> FeosResult<Pressure<Gradient<P>>>
224224
where
225225
Self: Residual<U2, Gradient<P>>,
226226
{
227227
let eos_f64 = self.re();
228+
let (liquid_molefracs, _) = liquid_molefracs.into_molefracs(&eos_f64)?;
228229
let vle = PhaseEquilibrium::bubble_point(
229230
&eos_f64,
230231
temperature,
231-
&liquid_molefracs,
232+
liquid_molefracs,
232233
pressure,
233234
None,
234235
Default::default(),
@@ -265,20 +266,21 @@ pub trait PropertiesAD {
265266
Ok(Pressure::from_reduced(p))
266267
}
267268

268-
fn dew_point_pressure<const P: usize>(
269+
fn dew_point_pressure<const P: usize, X: Composition<f64, U2>>(
269270
&self,
270271
temperature: Temperature,
271272
pressure: Option<Pressure>,
272-
vapor_molefracs: SVector<f64, 2>,
273+
vapor_molefracs: X,
273274
) -> FeosResult<Pressure<Gradient<P>>>
274275
where
275276
Self: Residual<U2, Gradient<P>>,
276277
{
277278
let eos_f64 = self.re();
279+
let (vapor_molefracs, _) = vapor_molefracs.into_molefracs(&eos_f64)?;
278280
let vle = PhaseEquilibrium::dew_point(
279281
&eos_f64,
280282
temperature,
281-
&vapor_molefracs,
283+
vapor_molefracs,
282284
pressure,
283285
None,
284286
Default::default(),
@@ -329,12 +331,8 @@ pub trait PropertiesAD {
329331
parameters,
330332
input,
331333
|eos: &Self::Lifted<Gradient<P>>, inp| {
332-
eos.bubble_point_pressure(
333-
inp[0] * KELVIN,
334-
Some(inp[2] * PASCAL),
335-
SVector::from([inp[1], 1.0 - inp[1]]),
336-
)
337-
.map(|p| p.convert_into(PASCAL))
334+
eos.bubble_point_pressure(inp[0] * KELVIN, Some(inp[2] * PASCAL), inp[1])
335+
.map(|p| p.convert_into(PASCAL))
338336
},
339337
)
340338
}
@@ -353,12 +351,8 @@ pub trait PropertiesAD {
353351
parameters,
354352
input,
355353
|eos: &Self::Lifted<Gradient<P>>, inp| {
356-
eos.dew_point_pressure(
357-
inp[0] * KELVIN,
358-
Some(inp[2] * PASCAL),
359-
SVector::from([inp[1], 1.0 - inp[1]]),
360-
)
361-
.map(|p| p.convert_into(PASCAL))
354+
eos.dew_point_pressure(inp[0] * KELVIN, Some(inp[2] * PASCAL), inp[1])
355+
.map(|p| p.convert_into(PASCAL))
362356
},
363357
)
364358
}

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

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,10 @@ pub enum FeosError {
2424
InvalidState(String, String, f64),
2525
#[error("Undetermined state: {0}")]
2626
UndeterminedState(String),
27+
#[error(
28+
"Extensive properties can only be evaluated for states that are initialized with extensive properties."
29+
)]
30+
IntensiveState,
2731
#[error("System is supercritical.")]
2832
SuperCritical,
2933
#[error("No phase split according to stability analysis.")]

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

Lines changed: 11 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -299,8 +299,8 @@ mod tests {
299299

300300
// residual properties
301301
assert_relative_eq!(
302-
s.helmholtz_energy(Contributions::Residual),
303-
sr.residual_helmholtz_energy(),
302+
s.helmholtz_energy(Contributions::Residual)?,
303+
sr.residual_helmholtz_energy()?,
304304
max_relative = 1e-15
305305
);
306306
assert_relative_eq!(
@@ -309,8 +309,8 @@ mod tests {
309309
max_relative = 1e-15
310310
);
311311
assert_relative_eq!(
312-
s.entropy(Contributions::Residual),
313-
sr.residual_entropy(),
312+
s.entropy(Contributions::Residual)?,
313+
sr.residual_entropy()?,
314314
max_relative = 1e-15
315315
);
316316
assert_relative_eq!(
@@ -319,8 +319,8 @@ mod tests {
319319
max_relative = 1e-15
320320
);
321321
assert_relative_eq!(
322-
s.enthalpy(Contributions::Residual),
323-
sr.residual_enthalpy(),
322+
s.enthalpy(Contributions::Residual)?,
323+
sr.residual_enthalpy()?,
324324
max_relative = 1e-15
325325
);
326326
assert_relative_eq!(
@@ -329,8 +329,8 @@ mod tests {
329329
max_relative = 1e-15
330330
);
331331
assert_relative_eq!(
332-
s.internal_energy(Contributions::Residual),
333-
sr.residual_internal_energy(),
332+
s.internal_energy(Contributions::Residual)?,
333+
sr.residual_internal_energy()?,
334334
max_relative = 1e-15
335335
);
336336
assert_relative_eq!(
@@ -339,12 +339,12 @@ mod tests {
339339
max_relative = 1e-15
340340
);
341341
assert_relative_eq!(
342-
s.gibbs_energy(Contributions::Residual)
343-
- s.total_moles()
342+
s.gibbs_energy(Contributions::Residual)?
343+
- s.total_moles()?
344344
* RGAS
345345
* s.temperature
346346
* s.compressibility(Contributions::Total).ln(),
347-
sr.residual_gibbs_energy(),
347+
sr.residual_gibbs_energy()?,
348348
max_relative = 1e-15
349349
);
350350
assert_relative_eq!(

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

Lines changed: 36 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,7 @@ use crate::state::{
44
Contributions,
55
DensityInitialization::{InitialDensity, Liquid, Vapor},
66
};
7-
use crate::{ReferenceSystem, Residual, SolverOptions, State, Verbosity};
7+
use crate::{Composition, ReferenceSystem, Residual, SolverOptions, State, Verbosity};
88
use nalgebra::allocator::Allocator;
99
use nalgebra::{DMatrix, DVector, DefaultAllocator, Dim, Dyn, OVector, U1};
1010
#[cfg(feature = "ndarray")]
@@ -141,10 +141,10 @@ where
141141
{
142142
/// Calculate a phase equilibrium for a given temperature
143143
/// or pressure and composition of the liquid phase.
144-
pub fn bubble_point<TP: TemperatureOrPressure<D>>(
144+
pub fn bubble_point<TP: TemperatureOrPressure<D>, X: Composition<D, N>>(
145145
eos: &E,
146146
temperature_or_pressure: TP,
147-
liquid_molefracs: &OVector<D, N>,
147+
liquid_molefracs: X,
148148
tp_init: Option<TP::Other>,
149149
vapor_molefracs: Option<&OVector<f64, N>>,
150150
options: (SolverOptions, SolverOptions),
@@ -162,10 +162,10 @@ where
162162

163163
/// Calculate a phase equilibrium for a given temperature
164164
/// or pressure and composition of the vapor phase.
165-
pub fn dew_point<TP: TemperatureOrPressure<D>>(
165+
pub fn dew_point<TP: TemperatureOrPressure<D>, X: Composition<D, N>>(
166166
eos: &E,
167167
temperature_or_pressure: TP,
168-
vapor_molefracs: &OVector<D, N>,
168+
vapor_molefracs: X,
169169
tp_init: Option<TP::Other>,
170170
liquid_molefracs: Option<&OVector<f64, N>>,
171171
options: (SolverOptions, SolverOptions),
@@ -181,35 +181,43 @@ where
181181
)
182182
}
183183

184-
pub(super) fn bubble_dew_point<TP: TemperatureOrPressure<D>>(
184+
pub(super) fn bubble_dew_point<TP: TemperatureOrPressure<D>, X: Composition<D, N>>(
185185
eos: &E,
186186
temperature_or_pressure: TP,
187-
vapor_molefracs: &OVector<D, N>,
187+
vapor_molefracs: X,
188188
tp_init: Option<TP::Other>,
189189
liquid_molefracs: Option<&OVector<f64, N>>,
190190
bubble: bool,
191191
options: (SolverOptions, SolverOptions),
192192
) -> FeosResult<Self> {
193-
let (temperature, pressure, iterate_p) =
194-
temperature_or_pressure.temperature_pressure(tp_init);
195-
Self::bubble_dew_point_tp(
196-
eos,
197-
temperature,
198-
pressure,
199-
vapor_molefracs,
200-
liquid_molefracs,
201-
bubble,
202-
iterate_p,
203-
options,
204-
)
193+
if eos.components() == 1 {
194+
let mut vle = Self::pure(eos, temperature_or_pressure, None, options.1)?;
195+
if bubble {
196+
vle.phase_fractions = [D::from(0.0), D::from(1.0)];
197+
}
198+
Ok(vle)
199+
} else {
200+
let (temperature, pressure, iterate_p) =
201+
temperature_or_pressure.temperature_pressure(tp_init);
202+
Self::bubble_dew_point_tp(
203+
eos,
204+
temperature,
205+
pressure,
206+
vapor_molefracs,
207+
liquid_molefracs,
208+
bubble,
209+
iterate_p,
210+
options,
211+
)
212+
}
205213
}
206214

207215
#[expect(clippy::too_many_arguments)]
208-
fn bubble_dew_point_tp(
216+
fn bubble_dew_point_tp<X: Composition<D, N>>(
209217
eos: &E,
210218
temperature: Option<Temperature<D>>,
211219
pressure: Option<Pressure<D>>,
212-
molefracs_spec: &OVector<D, N>,
220+
composition: X,
213221
molefracs_init: Option<&OVector<f64, N>>,
214222
bubble: bool,
215223
iterate_p: bool,
@@ -218,6 +226,7 @@ where
218226
let eos_re = eos.re();
219227
let mut temperature_re = temperature.map(|t| t.re());
220228
let mut pressure_re = pressure.map(|p| p.re());
229+
let (molefracs_spec, total_moles) = composition.into_molefracs(eos)?;
221230
let molefracs_spec_re = molefracs_spec.map(|x| x.re());
222231
let (v1, rho2) = if iterate_p {
223232
// temperature is specified
@@ -306,7 +315,7 @@ where
306315
Self::newton_step_t(
307316
eos,
308317
t,
309-
molefracs_spec,
318+
&molefracs_spec,
310319
&mut p,
311320
&mut molar_volume,
312321
&mut rho2,
@@ -316,7 +325,7 @@ where
316325
Self::newton_step_p(
317326
eos,
318327
&mut t,
319-
molefracs_spec,
328+
&molefracs_spec,
320329
p,
321330
&mut molar_volume,
322331
&mut rho2,
@@ -339,11 +348,11 @@ where
339348
x2,
340349
)?;
341350

342-
Ok(PhaseEquilibrium(if bubble {
343-
[state2, state1]
351+
Ok(if bubble {
352+
PhaseEquilibrium::with_vapor_phase_fraction(state2, state1, D::from(0.0), total_moles)
344353
} else {
345-
[state1, state2]
346-
}))
354+
PhaseEquilibrium::with_vapor_phase_fraction(state1, state2, D::from(1.0), total_moles)
355+
})
347356
}
348357

349358
fn newton_step_t(

0 commit comments

Comments
 (0)