Skip to content

Commit 1c41ed5

Browse files
prehnerprehner
committed
Automatically detect VLLEs in BinaryPhaseDiagram
1 parent d7d99fd commit 1c41ed5

File tree

1 file changed

+174
-7
lines changed

1 file changed

+174
-7
lines changed

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

Lines changed: 174 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,7 @@ use super::{PhaseDiagram, PhaseEquilibrium};
33
use crate::errors::{FeosError, FeosResult};
44
use crate::state::{Contributions, DensityInitialization::Vapor, State, StateBuilder};
55
use crate::{ReferenceSystem, Residual, SolverOptions, Subset};
6-
use nalgebra::{DVector, dvector, stack, vector};
6+
use nalgebra::{DVector, dvector, matrix, stack, vector};
77
use ndarray::{Array1, s};
88
use num_dual::linalg::LU;
99
use quantity::{Density, Moles, Pressure, RGAS, Temperature};
@@ -32,7 +32,7 @@ impl<E: Residual + Subset> PhaseDiagram<E, 2> {
3232

3333
// Only calculate up to specified compositions
3434
if let Some(x_lle) = x_lle {
35-
let (states1, states2) = Self::calculate_vlle(
35+
let [states1, states2] = Self::calculate_vlle(
3636
eos,
3737
temperature_or_pressure,
3838
npoints,
@@ -94,18 +94,18 @@ impl<E: Residual + Subset> PhaseDiagram<E, 2> {
9494
if !bubble {
9595
states = states.into_iter().rev().collect();
9696
}
97+
let states = check_for_vlle(temperature_or_pressure, states, npoints, bubble_dew_options);
9798
Ok(Self { states })
9899
}
99100

100-
#[expect(clippy::type_complexity)]
101101
fn calculate_vlle<TP: TemperatureOrPressure>(
102102
eos: &E,
103103
tp: TP,
104104
npoints: usize,
105105
x_lle: (f64, f64),
106106
vle_sat: [Option<PhaseEquilibrium<E, 2>>; 2],
107107
bubble_dew_options: (SolverOptions, SolverOptions),
108-
) -> FeosResult<(Vec<PhaseEquilibrium<E, 2>>, Vec<PhaseEquilibrium<E, 2>>)> {
108+
) -> FeosResult<[Vec<PhaseEquilibrium<E, 2>>; 2]> {
109109
match vle_sat {
110110
[Some(vle2), Some(vle1)] => {
111111
let states1 = iterate_vle(
@@ -128,7 +128,7 @@ impl<E: Residual + Subset> PhaseDiagram<E, 2> {
128128
true,
129129
bubble_dew_options,
130130
);
131-
Ok((states1, states2))
131+
Ok([states1, states2])
132132
}
133133
_ => Err(FeosError::SuperCritical),
134134
}
@@ -228,6 +228,174 @@ fn iterate_vle<E: Residual + Subset, TP: TemperatureOrPressure>(
228228
vle_vec
229229
}
230230

231+
fn check_for_vlle<E: Residual + Subset, TP: TemperatureOrPressure>(
232+
tp: TP,
233+
states: Vec<PhaseEquilibrium<E, 2>>,
234+
npoints: usize,
235+
bubble_dew_options: (SolverOptions, SolverOptions),
236+
) -> Vec<PhaseEquilibrium<E, 2>> {
237+
let n = states.len();
238+
let p: Vec<_> = states
239+
.iter()
240+
.map(|s| s.vapor().pressure(Contributions::Total))
241+
.collect();
242+
let t: Vec<_> = states.iter().map(|s| s.vapor().temperature).collect();
243+
let x: Vec<_> = states.iter().map(|s| s.liquid().molefracs[0]).collect();
244+
let y: Vec<_> = states.iter().map(|s| s.vapor().molefracs[0]).collect();
245+
246+
// Determine if the dew line intersects with itself
247+
if let Some(t) = tp.temperature()
248+
&& p[1] > p[0]
249+
&& p[n - 2] > p[n - 1]
250+
{
251+
let [mut i, mut j] = [0, n - 1];
252+
while i != j {
253+
if p[i] > p[j] {
254+
j -= 1;
255+
} else {
256+
i += 1
257+
}
258+
if y[j] < y[i] {
259+
// intersection found!
260+
let (xj, yj, pj) = if j == n - 2 {
261+
// Use Henry constant of component 2
262+
let k_inf = (states[n - 1].liquid().ln_phi() - states[n - 1].vapor().ln_phi())
263+
.map(f64::exp)[1];
264+
(
265+
[1.0, 1.0 - 1.0 / k_inf],
266+
[1.0, 0.0],
267+
[p[n - 1], p[n - 1] * (2.0 - 1.0 / k_inf)],
268+
)
269+
} else {
270+
// or interpolate linearly
271+
([x[j + 1], x[j]], [y[j + 1], y[j]], [p[j + 1], p[j]])
272+
};
273+
let (xi, yi, pi) = if i == 1 {
274+
// Use Henry constant of component 1
275+
let k_inf =
276+
(states[0].liquid().ln_phi() - states[0].vapor().ln_phi()).map(f64::exp)[0];
277+
(
278+
[0.0, 1.0 / k_inf],
279+
[0.0, 1.0],
280+
[p[0], p[0] * (2.0 - 1.0 / k_inf)],
281+
)
282+
} else {
283+
// or interpolate linearly
284+
([x[i - 1], x[i]], [y[i - 1], y[i]], [p[i - 1], p[i]])
285+
};
286+
// calculate intersection
287+
let a = matrix![yi[1] - yi[0], yj[0] - yj[1];
288+
(pi[1] - pi[0]).into_reduced(), (pj[0] - pj[1]).into_reduced()];
289+
let b = vector![yj[0] - yi[0], (pj[0] - pi[0]).into_reduced()];
290+
let [[r, s]] = LU::new(a).unwrap().solve(&b).data.0;
291+
let (xi, xj, p) = (
292+
xi[0] + r * (xi[1] - xi[0]),
293+
xj[0] + s * (xj[1] - xj[0]),
294+
pi[0] + r * (pi[1] - pi[0]),
295+
);
296+
let Ok(vlle) = PhaseEquilibrium::heteroazeotrope(
297+
&states[0].liquid().eos,
298+
t,
299+
(xi, xj),
300+
Some(p),
301+
Default::default(),
302+
bubble_dew_options,
303+
) else {
304+
return states;
305+
};
306+
let x_hetero = (vlle.liquid1().molefracs[0], vlle.liquid2().molefracs[0]);
307+
return PhaseDiagram::binary_vle(
308+
&states[0].liquid().eos,
309+
tp,
310+
Some(npoints),
311+
Some(x_hetero),
312+
bubble_dew_options,
313+
)
314+
.map_or(states, |dia| dia.states);
315+
}
316+
}
317+
} else if let Some(p) = tp.pressure()
318+
&& t[1] < t[0]
319+
&& t[n - 2] < t[n - 1]
320+
{
321+
let [mut i, mut j] = [0, n - 1];
322+
while i != j {
323+
if t[i] < t[j] {
324+
j -= 1;
325+
} else {
326+
i += 1
327+
}
328+
if y[j] < y[i] {
329+
// intersection found!
330+
let (xj, yj, tj) = if j == n - 2 {
331+
// Use Henry constant of component 2
332+
let vle = &states[n - 1];
333+
let k_inf = (vle.liquid().ln_phi() - vle.vapor().ln_phi()).map(f64::exp)[1];
334+
let dh = vle.vapor().residual_molar_enthalpy()
335+
- vle.liquid().residual_molar_enthalpy();
336+
let dv = 1.0 / vle.vapor().density - 1.0 / vle.liquid().density;
337+
let pdv_dh = (p * dv).convert_into(dh);
338+
(
339+
[1.0, 1.0 - 1.0 / k_inf],
340+
[1.0, 0.0],
341+
[t[n - 1], t[n - 1] * (1.0 - (k_inf - 1.0) / k_inf * pdv_dh)],
342+
)
343+
} else {
344+
// or interpolate linearly
345+
([x[j + 1], x[j]], [y[j + 1], y[j]], [t[j + 1], t[j]])
346+
};
347+
let (xi, yi, ti) = if i == 1 {
348+
// Use Henry constant of component 1
349+
let vle = &states[0];
350+
let k_inf = (vle.liquid().ln_phi() - vle.vapor().ln_phi()).map(f64::exp)[0];
351+
let dh = vle.vapor().residual_molar_enthalpy()
352+
- vle.liquid().residual_molar_enthalpy();
353+
let dv = 1.0 / vle.vapor().density - 1.0 / vle.liquid().density;
354+
let pdv_dh = (p * dv).convert_into(dh);
355+
(
356+
[0.0, 1.0 / k_inf],
357+
[0.0, 1.0],
358+
[t[0], t[0] * (1.0 - (k_inf - 1.0) / k_inf * pdv_dh)],
359+
)
360+
} else {
361+
// or interpolate linearly
362+
([x[i - 1], x[i]], [y[i - 1], y[i]], [t[i - 1], t[i]])
363+
};
364+
// calculate intersection
365+
let a = matrix![yi[1] - yi[0], yj[0] - yj[1];
366+
(ti[1] - ti[0]).into_reduced(), (tj[0] - tj[1]).into_reduced()];
367+
let b = vector![yj[0] - yi[0], (tj[0] - ti[0]).into_reduced()];
368+
let [[r, s]] = LU::new(a).unwrap().solve(&b).data.0;
369+
let (xi, xj, t) = (
370+
xi[0] + r * (xi[1] - xi[0]),
371+
xj[0] + s * (xj[1] - xj[0]),
372+
ti[0] + r * (ti[1] - ti[0]),
373+
);
374+
let Ok(vlle) = PhaseEquilibrium::heteroazeotrope(
375+
&states[0].liquid().eos,
376+
p,
377+
(xi, xj),
378+
Some(t),
379+
Default::default(),
380+
bubble_dew_options,
381+
) else {
382+
return states;
383+
};
384+
let x_hetero = (vlle.liquid1().molefracs[0], vlle.liquid2().molefracs[0]);
385+
return PhaseDiagram::binary_vle(
386+
&states[0].liquid().eos,
387+
tp,
388+
Some(npoints),
389+
Some(x_hetero),
390+
bubble_dew_options,
391+
)
392+
.map_or(states, |dia| dia.states);
393+
}
394+
}
395+
}
396+
states
397+
}
398+
231399
/// Phase diagram (Txy or pxy) for a system with heteroazeotropic phase behavior.
232400
pub struct PhaseDiagramHetero<E> {
233401
pub vle1: PhaseDiagram<E, 2>,
@@ -270,7 +438,7 @@ impl<E: Residual + Subset> PhaseDiagram<E, 2> {
270438
let x_hetero = (vlle.liquid1().molefracs[0], vlle.liquid2().molefracs[0]);
271439

272440
// calculate vapor liquid equilibria
273-
let (dia1, dia2) = PhaseDiagram::calculate_vlle(
441+
let [dia1, dia2] = PhaseDiagram::calculate_vlle(
274442
eos,
275443
temperature_or_pressure,
276444
npoints_vle,
@@ -334,7 +502,6 @@ impl<E: Residual> PhaseEquilibrium<E, 3> {
334502
) -> FeosResult<Self> {
335503
let (temperature, pressure, iterate_p) =
336504
temperature_or_pressure.temperature_pressure(tp_init);
337-
// let tp_init = tp_init.map(|tp_init| temperature_or_pressure.temperature_pressure(tp_init));
338505
if iterate_p {
339506
PhaseEquilibrium::heteroazeotrope_t(
340507
eos,

0 commit comments

Comments
 (0)