Skip to content

Commit 46cef76

Browse files
prehnerprehner
authored
More robust phase equilibrium calculations (#23)
* More robust phase equilibrium calculations * minor refinements and updated CHANGELOG * add initial temperature/pressure as parameter for VLLE calculations * add updated example notebook
1 parent 052c291 commit 46cef76

File tree

12 files changed

+722
-189
lines changed

12 files changed

+722
-189
lines changed

‎examples/pcsaft/PhaseDiagramBinary.ipynb‎

Lines changed: 6 additions & 6 deletions
Large diffs are not rendered by default.

‎feos-core/CHANGELOG.md‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -5,9 +5,17 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
55
and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
66

77
## Unreleased
8+
### Added
9+
- Added `State::spinodal` that calculates both spinodal points for a given temperature and composition using the same objective function that is also used in the critical point calculation. [#23](https://github.com/feos-org/feos/pull/23)
10+
- Added `PhaseDiagram::bubble_point_line`, `PhaseDiagram::dew_point_line`, and `PhaseDiagram::spinodal` to calculate phase envelopes for mixtures with fixed compositions. [#23](https://github.com/feos-org/feos/pull/23)
11+
812
### Changed
913
- Made binary parameters in `from_records` Python routine an `Option`. [#35](https://github.com/feos-org/feos/pull/35)
1014
- Added panic with message when parsing missing Identifiers variants. [#35](https://github.com/feos-org/feos/pull/35)
15+
- Generalized the initialization routines for pure component VLEs at given temperature to multicomponent systems. [#23](https://github.com/feos-org/feos/pull/23)
16+
17+
### Removed
18+
- Removed the (internal) `SpinodalPoint` struct that was used within density iterations in favor of a simpler interface. [#23](https://github.com/feos-org/feos/pull/23)
1119

1220
### Fixed
1321
- Avoid panicking when calculating `ResidualNpt` properties of states with negative pressures. [#42](https://github.com/feos-org/feos/pull/42)

‎feos-core/src/density_iteration.rs‎

Lines changed: 20 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -5,12 +5,6 @@ use crate::EosUnit;
55
use quantity::{QuantityArray1, QuantityScalar};
66
use std::rc::Rc;
77

8-
pub struct SpinodalPoint<U: EosUnit> {
9-
pub p: QuantityScalar<U>,
10-
pub dp_drho: QuantityScalar<U>,
11-
pub rho: QuantityScalar<U>,
12-
}
13-
148
pub fn density_iteration<U: EosUnit, E: EquationOfState>(
159
eos: &Rc<E>,
1610
temperature: QuantityScalar<U>,
@@ -64,9 +58,9 @@ pub fn density_iteration<U: EosUnit, E: EquationOfState>(
6458
.2;
6559

6660
if rho > 0.85 * maxdensity {
67-
let sp = pressure_spinodal(eos, temperature, initial_density, moles)?;
68-
rho = sp.rho;
69-
error = sp.p - pressure;
61+
let [sp_p, sp_rho] = pressure_spinodal(eos, temperature, initial_density, moles)?;
62+
rho = sp_rho;
63+
error = sp_p - pressure;
7064
if rho > 0.85 * maxdensity {
7165
if error.is_sign_negative() {
7266
return Err(EosError::IterationFailed(String::from("density_iteration")));
@@ -79,29 +73,28 @@ pub fn density_iteration<U: EosUnit, E: EquationOfState>(
7973
rho = (rho * 1.1).min(maxdensity)?
8074
}
8175
} else if error.is_sign_positive() && d2pdrho2.is_sign_positive() {
82-
let sp = pressure_spinodal(eos, temperature, initial_density, moles)?;
83-
rho = sp.rho;
84-
error = sp.p - pressure;
76+
let [sp_p, sp_rho] = pressure_spinodal(eos, temperature, initial_density, moles)?;
77+
rho = sp_rho;
78+
error = sp_p - pressure;
8579
if error.is_sign_positive() {
8680
rho = 0.001 * maxdensity
8781
} else {
8882
rho = (rho * 1.1).min(maxdensity)?
8983
}
9084
} else if error.is_sign_negative() && d2pdrho2.is_sign_negative() {
91-
let sp = pressure_spinodal(eos, temperature, initial_density, moles)?;
92-
rho = sp.rho;
93-
error = sp.p - pressure;
85+
let [sp_p, sp_rho] = pressure_spinodal(eos, temperature, initial_density, moles)?;
86+
rho = sp_rho;
87+
error = sp_p - pressure;
9488
if error.is_sign_negative() {
9589
rho = 0.8 * maxdensity
9690
} else {
9791
rho = rho * 0.8
9892
}
9993
} else if error.is_sign_negative() && d2pdrho2.is_sign_positive() {
100-
let sp_l = pressure_spinodal(eos, temperature, 0.8 * maxdensity, moles)?;
101-
let rho_l = sp_l.rho;
102-
let sp_v = pressure_spinodal(eos, temperature, 0.001 * maxdensity, moles)?;
103-
let rho_v = sp_v.rho;
104-
error = sp_v.p - pressure;
94+
let [_, rho_l] = pressure_spinodal(eos, temperature, 0.8 * maxdensity, moles)?;
95+
let [sp_v_p, rho_v] =
96+
pressure_spinodal(eos, temperature, 0.001 * maxdensity, moles)?;
97+
error = sp_v_p - pressure;
10598
if error.is_sign_positive()
10699
&& (initial_density - rho_v).abs() < (initial_density - rho_l).abs()
107100
{
@@ -110,11 +103,10 @@ pub fn density_iteration<U: EosUnit, E: EquationOfState>(
110103
rho = (rho_l * 1.1).min(maxdensity)?
111104
}
112105
} else if error.is_sign_positive() && d2pdrho2.is_sign_negative() {
113-
let sp_l = pressure_spinodal(eos, temperature, 0.8 * maxdensity, moles)?;
114-
let rho_l = sp_l.rho;
115-
let sp_v = pressure_spinodal(eos, temperature, 0.001 * maxdensity, moles)?;
116-
let rho_v = sp_v.rho;
117-
error = sp_v.p - pressure;
106+
let [_, rho_l] = pressure_spinodal(eos, temperature, 0.8 * maxdensity, moles)?;
107+
let [sp_v_p, rho_v] =
108+
pressure_spinodal(eos, temperature, 0.001 * maxdensity, moles)?;
109+
error = sp_v_p - pressure;
118110
if error.is_sign_negative()
119111
&& (initial_density - rho_v).abs() > (initial_density - rho_l).abs()
120112
{
@@ -149,12 +141,12 @@ pub fn density_iteration<U: EosUnit, E: EquationOfState>(
149141
}
150142
}
151143

152-
pub fn pressure_spinodal<U: EosUnit, E: EquationOfState>(
144+
fn pressure_spinodal<U: EosUnit, E: EquationOfState>(
153145
eos: &Rc<E>,
154146
temperature: QuantityScalar<U>,
155147
rho_init: QuantityScalar<U>,
156148
moles: &QuantityArray1<U>,
157-
) -> EosResult<SpinodalPoint<U>> {
149+
) -> EosResult<[QuantityScalar<U>; 2]> {
158150
let maxiter = 30;
159151
let abstol = 1e-8;
160152

@@ -186,11 +178,7 @@ pub fn pressure_spinodal<U: EosUnit, E: EquationOfState>(
186178
.abs()
187179
< abstol
188180
{
189-
return Ok(SpinodalPoint {
190-
p,
191-
dp_drho: dpdrho,
192-
rho,
193-
});
181+
return Ok([p, rho]);
194182
}
195183
}
196184
Err(EosError::NotConverged("pressure_spinodal".to_owned()))

0 commit comments

Comments
 (0)