Skip to content

More robust phase equilibrium calculations #23

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
prehner merged 5 commits into from
Aug 25, 2022
Merged

More robust phase equilibrium calculations #23

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
0047a68
More robust phase equilibrium calculations
prehner Jun 17, 2022
1ef5412
minor refinements and updated CHANGELOG
prehner Aug 8, 2022
20ee2a8
add initial temperature/pressure as parameter for VLLE calculations
prehner Aug 8, 2022
acec37f
add updated example notebook
prehner Aug 8, 2022
5f2f039
Merge branch 'main' into robust_vle
prehner Aug 8, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 6 additions & 6 deletions examples/pcsaft/PhaseDiagramBinary.ipynb
View file
Open in desktop

Large diffs are not rendered by default.

8 changes: 8 additions & 0 deletions feos-core/CHANGELOG.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,9 +5,17 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).

## Unreleased
### Added
- 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)
- 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)

### Changed
- Made binary parameters in `from_records` Python routine an `Option`. [#35](https://github.com/feos-org/feos/pull/35)
- Added panic with message when parsing missing Identifiers variants. [#35](https://github.com/feos-org/feos/pull/35)
- Generalized the initialization routines for pure component VLEs at given temperature to multicomponent systems. [#23](https://github.com/feos-org/feos/pull/23)

### Removed
- 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)

### Fixed
- Avoid panicking when calculating `ResidualNpt` properties of states with negative pressures. [#42](https://github.com/feos-org/feos/pull/42)
Expand Down
52 changes: 20 additions & 32 deletions feos-core/src/density_iteration.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,12 +5,6 @@ use crate::EosUnit;
use quantity::{QuantityArray1, QuantityScalar};
use std::rc::Rc;

pub struct SpinodalPoint<U: EosUnit> {
pub p: QuantityScalar<U>,
pub dp_drho: QuantityScalar<U>,
pub rho: QuantityScalar<U>,
}

pub fn density_iteration<U: EosUnit, E: EquationOfState>(
eos: &Rc<E>,
temperature: QuantityScalar<U>,
Expand Down Expand Up @@ -64,9 +58,9 @@ pub fn density_iteration<U: EosUnit, E: EquationOfState>(
.2;

if rho > 0.85 * maxdensity {
let sp = pressure_spinodal(eos, temperature, initial_density, moles)?;
rho = sp.rho;
error = sp.p - pressure;
let [sp_p, sp_rho] = pressure_spinodal(eos, temperature, initial_density, moles)?;
rho = sp_rho;
error = sp_p - pressure;
if rho > 0.85 * maxdensity {
if error.is_sign_negative() {
return Err(EosError::IterationFailed(String::from("density_iteration")));
Expand All @@ -79,29 +73,28 @@ pub fn density_iteration<U: EosUnit, E: EquationOfState>(
rho = (rho * 1.1).min(maxdensity)?
}
} else if error.is_sign_positive() && d2pdrho2.is_sign_positive() {
let sp = pressure_spinodal(eos, temperature, initial_density, moles)?;
rho = sp.rho;
error = sp.p - pressure;
let [sp_p, sp_rho] = pressure_spinodal(eos, temperature, initial_density, moles)?;
rho = sp_rho;
error = sp_p - pressure;
if error.is_sign_positive() {
rho = 0.001 * maxdensity
} else {
rho = (rho * 1.1).min(maxdensity)?
}
} else if error.is_sign_negative() && d2pdrho2.is_sign_negative() {
let sp = pressure_spinodal(eos, temperature, initial_density, moles)?;
rho = sp.rho;
error = sp.p - pressure;
let [sp_p, sp_rho] = pressure_spinodal(eos, temperature, initial_density, moles)?;
rho = sp_rho;
error = sp_p - pressure;
if error.is_sign_negative() {
rho = 0.8 * maxdensity
} else {
rho = rho * 0.8
}
} else if error.is_sign_negative() && d2pdrho2.is_sign_positive() {
let sp_l = pressure_spinodal(eos, temperature, 0.8 * maxdensity, moles)?;
let rho_l = sp_l.rho;
let sp_v = pressure_spinodal(eos, temperature, 0.001 * maxdensity, moles)?;
let rho_v = sp_v.rho;
error = sp_v.p - pressure;
let [_, rho_l] = pressure_spinodal(eos, temperature, 0.8 * maxdensity, moles)?;
let [sp_v_p, rho_v] =
pressure_spinodal(eos, temperature, 0.001 * maxdensity, moles)?;
error = sp_v_p - pressure;
if error.is_sign_positive()
&& (initial_density - rho_v).abs() < (initial_density - rho_l).abs()
{
Expand All @@ -110,11 +103,10 @@ pub fn density_iteration<U: EosUnit, E: EquationOfState>(
rho = (rho_l * 1.1).min(maxdensity)?
}
} else if error.is_sign_positive() && d2pdrho2.is_sign_negative() {
let sp_l = pressure_spinodal(eos, temperature, 0.8 * maxdensity, moles)?;
let rho_l = sp_l.rho;
let sp_v = pressure_spinodal(eos, temperature, 0.001 * maxdensity, moles)?;
let rho_v = sp_v.rho;
error = sp_v.p - pressure;
let [_, rho_l] = pressure_spinodal(eos, temperature, 0.8 * maxdensity, moles)?;
let [sp_v_p, rho_v] =
pressure_spinodal(eos, temperature, 0.001 * maxdensity, moles)?;
error = sp_v_p - pressure;
if error.is_sign_negative()
&& (initial_density - rho_v).abs() > (initial_density - rho_l).abs()
{
Expand Down Expand Up @@ -149,12 +141,12 @@ pub fn density_iteration<U: EosUnit, E: EquationOfState>(
}
}

pub fn pressure_spinodal<U: EosUnit, E: EquationOfState>(
fn pressure_spinodal<U: EosUnit, E: EquationOfState>(
eos: &Rc<E>,
temperature: QuantityScalar<U>,
rho_init: QuantityScalar<U>,
moles: &QuantityArray1<U>,
) -> EosResult<SpinodalPoint<U>> {
) -> EosResult<[QuantityScalar<U>; 2]> {
let maxiter = 30;
let abstol = 1e-8;

Expand Down Expand Up @@ -186,11 +178,7 @@ pub fn pressure_spinodal<U: EosUnit, E: EquationOfState>(
.abs()
< abstol
{
return Ok(SpinodalPoint {
p,
dp_drho: dpdrho,
rho,
});
return Ok([p, rho]);
}
}
Err(EosError::NotConverged("pressure_spinodal".to_owned()))
Expand Down
Loading