Skip to content

Commit 2b9f8e4

Browse files
prehnerprehner
committed
include phase fraction in PhaseEquilibrium
1 parent 29be838 commit 2b9f8e4

File tree

10 files changed

+161
-107
lines changed

10 files changed

+161
-107
lines changed

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

Lines changed: 26 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -190,18 +190,28 @@ where
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+
Self::pure_with(
195+
eos,
196+
temperature_or_pressure,
197+
D::from(if bubble { 0.0 } else { 1.0 }),
198+
None,
199+
options.1,
200+
)
201+
} else {
202+
let (temperature, pressure, iterate_p) =
203+
temperature_or_pressure.temperature_pressure(tp_init);
204+
Self::bubble_dew_point_tp(
205+
eos,
206+
temperature,
207+
pressure,
208+
vapor_molefracs,
209+
liquid_molefracs,
210+
bubble,
211+
iterate_p,
212+
options,
213+
)
214+
}
205215
}
206216

207217
#[expect(clippy::too_many_arguments)]
@@ -339,11 +349,11 @@ where
339349
x2,
340350
)?;
341351

342-
Ok(PhaseEquilibrium(if bubble {
343-
[state2, state1]
352+
Ok(if bubble {
353+
PhaseEquilibrium([state2, state1], [D::from(0.0), D::from(1.0)])
344354
} else {
345-
[state1, state2]
346-
}))
355+
PhaseEquilibrium([state1, state2], [D::from(1.0), D::from(0.0)])
356+
})
347357
}
348358

349359
fn newton_step_t(

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

Lines changed: 51 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -4,8 +4,8 @@ use crate::state::{DensityInitialization, State};
44
use crate::{Contributions, ReferenceSystem};
55
use nalgebra::allocator::Allocator;
66
use nalgebra::{DefaultAllocator, Dim, Dyn, OVector};
7-
use num_dual::{DualNum, DualStruct, Gradients};
8-
use quantity::{Energy, Moles, Pressure, RGAS, Temperature};
7+
use num_dual::{DualNum, Gradients};
8+
use quantity::{MolarEnergy, Pressure, RGAS, Temperature};
99
use std::fmt;
1010
use std::fmt::Write;
1111

@@ -40,6 +40,7 @@ pub use phase_diagram_pure::PhaseDiagram;
4040
#[derive(Debug, Clone)]
4141
pub struct PhaseEquilibrium<E, const P: usize, N: Dim = Dyn, D: DualNum<f64> + Copy = f64>(
4242
pub [State<E, N, D>; P],
43+
pub [D; P],
4344
)
4445
where
4546
DefaultAllocator: Allocator<N>;
@@ -113,19 +114,38 @@ impl<E> PhaseEquilibrium<E, 3> {
113114
}
114115
}
115116

116-
impl<E: Residual<N>, N: Dim> PhaseEquilibrium<E, 2, N>
117+
impl<E: Residual<N, D>, N: Dim, D: DualNum<f64> + Copy> PhaseEquilibrium<E, 2, N, D>
117118
where
118119
DefaultAllocator: Allocator<N>,
119120
{
120-
pub(super) fn from_states(state1: State<E, N>, state2: State<E, N>) -> Self {
121-
let (vapor, liquid) = if state1.density.re() < state2.density.re() {
122-
(state1, state2)
123-
} else {
124-
(state2, state1)
125-
};
126-
Self([vapor, liquid])
121+
pub(super) fn single_phase(state: &State<E, N, D>) -> Self {
122+
Self::two_phase(state.clone(), state.clone())
127123
}
128124

125+
pub(super) fn two_phase(vapor: State<E, N, D>, liquid: State<E, N, D>) -> Self {
126+
Self::with_vapor_phase_fraction(vapor, liquid, D::from(1.0))
127+
}
128+
129+
pub(super) fn with_vapor_phase_fraction(
130+
vapor: State<E, N, D>,
131+
liquid: State<E, N, D>,
132+
vapor_phase_fraction: D,
133+
) -> Self {
134+
Self(
135+
[vapor, liquid],
136+
[vapor_phase_fraction, -vapor_phase_fraction + 1.0],
137+
)
138+
}
139+
140+
// fn from_states(state1: State<E, N>, state2: State<E, N>, beta: f64) -> Self {
141+
// let (vapor, liquid) = if state1.density.re() < state2.density.re() {
142+
// (state1, state2)
143+
// } else {
144+
// (state2, state1)
145+
// };
146+
// Self([vapor, liquid], [beta, 1.0 - beta])
147+
// }
148+
129149
// /// Creates a new PhaseEquilibrium that contains two states at the
130150
// /// specified temperature, pressure and molefracs.
131151
// ///
@@ -156,12 +176,12 @@ where
156176
// Ok(Self([vapor, liquid]))
157177
// }
158178

159-
pub(super) fn vapor_phase_fraction(&self) -> Option<f64> {
160-
self.vapor()
161-
.total_moles
162-
.zip(self.liquid().total_moles)
163-
.map(|(v, l)| (v / (l + v)).into_value())
164-
}
179+
// pub(super) fn vapor_phase_fraction(&self) -> Option<f64> {
180+
// self.vapor()
181+
// .total_moles
182+
// .zip(self.liquid().total_moles)
183+
// .map(|(v, l)| (v / (l + v)).into_value())
184+
// }
165185
}
166186

167187
impl<E: Residual<N>, N: Gradients, const P: usize> PhaseEquilibrium<E, P, N>
@@ -185,31 +205,35 @@ where
185205
Ok(self)
186206
}
187207

188-
pub(super) fn update_moles(
208+
pub(super) fn update_composition(
189209
&mut self,
190210
pressure: Pressure,
191-
moles: [&Moles<OVector<f64, N>>; P],
211+
molefracs: [&OVector<f64, N>; P],
212+
phase_fractions: [f64; P],
192213
) -> FeosResult<()> {
193214
for (i, s) in self.0.iter_mut().enumerate() {
194215
*s = State::new_npt(
195216
&s.eos,
196217
s.temperature,
197218
pressure,
198-
moles[i],
219+
molefracs[i],
199220
Some(DensityInitialization::InitialDensity(s.density)),
200221
)?;
201222
}
223+
self.1 = phase_fractions;
202224
Ok(())
203225
}
204226

205-
// Total Gibbs energy excluding the constant contribution RT sum_i N_i ln(\Lambda_i^3)
206-
pub(super) fn total_gibbs_energy(&self) -> Energy {
207-
self.0.iter().fold(Energy::from_reduced(0.0), |acc, s| {
208-
let ln_rho_m1 = s.partial_density().to_reduced().map(|r| r.ln() - 1.0);
209-
acc + s.residual_helmholtz_energy()
210-
+ s.pressure(Contributions::Total) * s.volume()
211-
+ RGAS * s.temperature * s.total_moles() * s.molefracs.dot(&ln_rho_m1)
212-
})
227+
// Total molar Gibbs energy excluding the constant contribution RT sum_i x_i ln(\Lambda_i^3)
228+
pub(super) fn total_molar_gibbs_energy(&self) -> MolarEnergy {
229+
self.0
230+
.iter()
231+
.fold(MolarEnergy::from_reduced(0.0), |acc, s| {
232+
let ln_rho_m1 = s.partial_density().to_reduced().map(|r| r.ln() - 1.0);
233+
acc + s.residual_molar_helmholtz_energy()
234+
+ s.pressure(Contributions::Total) * s.molar_volume
235+
+ RGAS * s.temperature * s.molefracs.dot(&ln_rho_m1)
236+
})
213237
}
214238
}
215239

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

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -63,7 +63,7 @@ impl<E: Residual + Subset> PhaseDiagram<E, 2> {
6363
None,
6464
SolverOptions::default(),
6565
)?;
66-
let cp_vle = PhaseEquilibrium::from_states(cp.clone(), cp.clone());
66+
let cp_vle = PhaseEquilibrium::single_phase(&cp);
6767
([0.0, cp.molefracs[0]], (vle2, cp_vle), bubble)
6868
}
6969
[None, Some(vle1)] => {
@@ -75,7 +75,7 @@ impl<E: Residual + Subset> PhaseDiagram<E, 2> {
7575
None,
7676
SolverOptions::default(),
7777
)?;
78-
let cp_vle = PhaseEquilibrium::from_states(cp.clone(), cp.clone());
78+
let cp_vle = PhaseEquilibrium::single_phase(&cp);
7979
([1.0, cp.molefracs[0]], (vle1, cp_vle), bubble)
8080
}
8181
[Some(vle2), Some(vle1)] => ([0.0, 1.0], (vle2, vle1), true),
@@ -437,7 +437,7 @@ impl<E: Residual> PhaseEquilibrium<E, 3> {
437437

438438
// check for convergence
439439
if res.norm() < options.tol.unwrap_or(TOL_HETERO) {
440-
return Ok(Self([v, l1, l2]));
440+
return Ok(Self([v, l1, l2], [1.0, 0.0, 0.0]));
441441
}
442442

443443
// calculate Jacobian
@@ -557,7 +557,7 @@ impl<E: Residual> PhaseEquilibrium<E, 3> {
557557

558558
// check for convergence
559559
if res.norm() < options.tol.unwrap_or(TOL_HETERO) {
560-
return Ok(Self([v, l1, l2]));
560+
return Ok(Self([v, l1, l2], [1.0, 0.0, 0.0]));
561561
}
562562

563563
let jacobian = stack![

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

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -54,7 +54,7 @@ impl<E: Residual> PhaseDiagram<E, 2> {
5454
states.push(vle.clone());
5555
}
5656
}
57-
states.push(PhaseEquilibrium::from_states(sc.clone(), sc));
57+
states.push(PhaseEquilibrium::single_phase(&sc));
5858

5959
Ok(PhaseDiagram::new(states))
6060
}

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

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -51,7 +51,7 @@ impl<E: Residual> PhaseDiagram<E, 2> {
5151
states.push(vle.clone());
5252
}
5353
}
54-
states.push(PhaseEquilibrium::from_states(sc.clone(), sc));
54+
states.push(PhaseEquilibrium::single_phase(&sc));
5555

5656
Ok(PhaseDiagram::new(states))
5757
}
@@ -116,7 +116,7 @@ impl<E: Residual> PhaseDiagram<E, 2> {
116116
}
117117
}
118118

119-
states.push(PhaseEquilibrium::from_states(sc.clone(), sc));
119+
states.push(PhaseEquilibrium::single_phase(&sc));
120120

121121
Ok(PhaseDiagram::new(states))
122122
}
@@ -146,11 +146,11 @@ impl<E: Residual> PhaseDiagram<E, 2> {
146146

147147
for ti in &temperatures {
148148
let spinodal = State::spinodal(eos, ti, molefracs, options).ok();
149-
if let Some(spinodal) = spinodal {
150-
states.push(PhaseEquilibrium(spinodal));
149+
if let Some([sp_v, sp_l]) = spinodal {
150+
states.push(PhaseEquilibrium::two_phase(sp_v, sp_l));
151151
}
152152
}
153-
states.push(PhaseEquilibrium::from_states(sc.clone(), sc));
153+
states.push(PhaseEquilibrium::single_phase(&sc));
154154

155155
Ok(PhaseDiagram::new(states))
156156
}

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

Lines changed: 26 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@ use nalgebra::{DefaultAllocator, Dim, Matrix3, OVector, SVector, U1, U2, vector}
88
use num_dual::{
99
Dual, Dual2Vec, DualNum, DualStruct, Gradients, first_derivative, implicit_derivative_sp,
1010
};
11-
use quantity::{Dimensionless, MOL, MolarVolume, Moles, Pressure, Quantity, Temperature};
11+
use quantity::{MOL, MolarVolume, Moles, Pressure, Quantity, Temperature};
1212

1313
const MAX_ITER_TP: usize = 400;
1414
const TOL_TP: f64 = 1e-8;
@@ -54,15 +54,14 @@ impl<E: Residual<U2, D>, D: DualNum<f64> + Copy> PhaseEquilibrium<E, 2, U2, D> {
5454
) -> FeosResult<Self> {
5555
let z = feed.get(0).convert_into(feed.get(0) + feed.get(1));
5656
let total_moles = feed.sum();
57-
let moles = vector![z.re(), 1.0 - z.re()] * MOL;
58-
let vle_re = State::new_npt(&eos.re(), temperature.re(), pressure.re(), moles, None)?
57+
let vle_re = State::new_npt(&eos.re(), temperature.re(), pressure.re(), z.re(), None)?
5958
.tp_flash(None, options, None)?;
6059

6160
// implicit differentiation
6261
let t = temperature.into_reduced();
6362
let p = pressure.into_reduced();
64-
let v_l = vle_re.liquid().density.into_reduced().recip();
65-
let v_v = vle_re.vapor().density.into_reduced().recip();
63+
let v_l = vle_re.liquid().molar_volume.into_reduced();
64+
let v_v = vle_re.vapor().molar_volume.into_reduced();
6665
let x = vle_re.liquid().molefracs[0];
6766
let y = vle_re.vapor().molefracs[0];
6867
let x = implicit_derivative_sp(
@@ -88,7 +87,10 @@ impl<E: Residual<U2, D>, D: DualNum<f64> + Copy> PhaseEquilibrium<E, 2, U2, D> {
8887
let moles = Quantity::new(vector![x, -x + 1.0] * phi * total_moles.convert_into(MOL));
8988
State::new_nvt(eos, temperature, volume, moles)
9089
};
91-
Ok(Self([state(y, v_v, beta)?, state(x, v_l, -beta + 1.0)?]))
90+
Ok(Self(
91+
[state(y, v_v, beta)?, state(x, v_l, -beta + 1.0)?],
92+
[beta, -beta + 1.0],
93+
))
9294
}
9395
}
9496

@@ -173,12 +175,12 @@ where
173175
// check convergence
174176
// unwrap is safe here, because after the first successive substitution step the
175177
// phase amounts in new_vle_state are known.
176-
let beta = new_vle_state.vapor_phase_fraction().unwrap();
177178
let tpd = [
178179
self.tangent_plane_distance(new_vle_state.vapor()),
179180
self.tangent_plane_distance(new_vle_state.liquid()),
180181
];
181-
let dg = (1.0 - beta) * tpd[1] + beta * tpd[0];
182+
let b = new_vle_state.1;
183+
let dg = b[0] * tpd[0] + b[1] * tpd[1];
182184

183185
// fix if only tpd[1] is positive
184186
if tpd[0] < 0.0 && dg >= 0.0 {
@@ -276,7 +278,7 @@ where
276278
}
277279

278280
// calculate total Gibbs energy before the extrapolation
279-
let gibbs = self.total_gibbs_energy();
281+
let gibbs = self.total_molar_gibbs_energy();
280282

281283
// extrapolate K values
282284
let delta_vec = [
@@ -304,7 +306,7 @@ where
304306
// calculate new states
305307
let mut trial_vle_state = self.clone();
306308
trial_vle_state.update_states(feed_state, &k)?;
307-
if trial_vle_state.total_gibbs_energy() < gibbs {
309+
if trial_vle_state.total_molar_gibbs_energy() < gibbs {
308310
*self = trial_vle_state;
309311
}
310312
}
@@ -366,23 +368,21 @@ where
366368

367369
fn update_states(&mut self, feed_state: &State<E, N>, k: &OVector<f64, N>) -> FeosResult<()> {
368370
// calculate vapor phase fraction using Rachford-Rice algorithm
369-
let beta = self.vapor_phase_fraction();
370-
let beta = rachford_rice(&feed_state.molefracs, k, beta)?;
371+
let beta = self.1[0];
372+
let b = rachford_rice(&feed_state.molefracs, k, Some(beta))?;
371373

372374
// update VLE
373375
let v = feed_state
374-
.moles()
375-
.clone()
376-
.component_mul(&Dimensionless::new(
377-
k.map(|k| beta * k / (1.0 - beta + beta * k)),
378-
));
376+
.molefracs
377+
.component_mul(&k.map(|k| b * k / (1.0 - b + b * k)));
379378
let l = feed_state
380-
.moles()
381-
.clone()
382-
.component_mul(&Dimensionless::new(
383-
k.map(|k| (1.0 - beta) / (1.0 - beta + beta * k)),
384-
));
385-
self.update_moles(feed_state.pressure(Contributions::Total), [&v, &l])?;
379+
.molefracs
380+
.component_mul(&k.map(|k| (1.0 - b) / (1.0 - b + b * k)));
381+
self.update_composition(
382+
feed_state.pressure(Contributions::Total),
383+
[&v, &l],
384+
[b, 1.0 - b],
385+
)?;
386386
Ok(())
387387
}
388388

@@ -391,9 +391,10 @@ where
391391
let state1 = stable_states.pop();
392392
let state2 = stable_states.pop();
393393
if let Some(s1) = state1 {
394-
let init1 = Self::from_states(s1.clone(), feed_state.clone());
394+
// TODO: sort
395+
let init1 = Self::two_phase(s1.clone(), feed_state.clone());
395396
if let Some(s2) = state2 {
396-
Ok((Self::from_states(s1, s2), Some(init1)))
397+
Ok((Self::two_phase(s1, s2), Some(init1)))
397398
} else {
398399
Ok((init1, None))
399400
}

0 commit comments

Comments
 (0)