From 268136fa0235a79e3b3a726b405ead28bb2fd1b8 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Sun, 15 Mar 2026 14:12:22 +0100 Subject: [PATCH] Implement ph and ps flashes for binary mixtures --- CHANGELOG.md | 5 +- crates/feos-core/src/equation_of_state/mod.rs | 115 +++-- crates/feos-core/src/lib.rs | 4 +- crates/feos-core/src/phase_equilibria/mod.rs | 42 +- .../src/phase_equilibria/px_flashes.rs | 465 ++++++++++++++++++ crates/feos-core/src/state/composition.rs | 31 +- crates/feos-core/src/state/mod.rs | 66 +-- crates/feos-core/src/state/properties.rs | 18 +- crates/feos-derive/src/ideal_gas.rs | 2 +- crates/feos-dft/src/profile/properties.rs | 2 +- crates/feos/src/ideal_gas/joback.rs | 34 +- crates/feos/src/multiparameter/mod.rs | 4 +- crates/feos/tests/pcsaft/mod.rs | 1 + crates/feos/tests/pcsaft/px_flashes.rs | 138 ++++++ crates/feos/tests/pcsaft/tp_flash.rs | 14 +- docs/recipes/index.md | 1 + .../recipes_phase_equilibrium_flash.ipynb | 149 ++++++ py-feos/src/phase_equilibria.rs | 152 ++++++ 18 files changed, 1110 insertions(+), 133 deletions(-) create mode 100644 crates/feos-core/src/phase_equilibria/px_flashes.rs create mode 100644 crates/feos/tests/pcsaft/px_flashes.rs create mode 100644 docs/recipes/recipes_phase_equilibrium_flash.ipynb diff --git a/CHANGELOG.md b/CHANGELOG.md index 513a7cc1e..3d69100a7 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,15 +10,18 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Rewrote `PhaseEquilibrium::pure_p` to mirror `pure_t` and enabled automatic differentiation. [#337](https://github.com/feos-org/feos/pull/337) - Added `boiling_temperature` to the list of properties for parallel evaluations of gradients. [#337](https://github.com/feos-org/feos/pull/337) - Added the `Composition` trait to allow more flexibility in the creation of states and phase equilibria. [#330](https://github.com/feos-org/feos/pull/330) +- Added `PhaseEquilibrium::ph_flash` and `PhaseEquilibrium::ps_flash`. [#338](https://github.com/feos-org/feos/pull/338) +- Added getters for `vapor_phase_fraction`, `molar_enthalpy`, `molar_entropy`, `total_moles`, `enthalpy`, and `entropy` to `PhaseEquilibrium`. [#338](https://github.com/feos-org/feos/pull/338) ### Changed - Removed any assumptions about the total number of moles in a `State` or `PhaseEquilibrium`. Evaluating extensive properties now returns a `Result`. [#330](https://github.com/feos-org/feos/pull/330) +- Redesigned the `IdealGas` trait and added `IdealGasAD` in analogy to `ResidualDyn` and `Residual`. [#330](https://github.com/feos-org/feos/pull/330) ### Removed - Removed the `StateBuilder` struct, because it is mostly obsolete with the addition of the `Composition` trait. [#330](https://github.com/feos-org/feos/pull/330) ### Packaging -- Updated `quantity` dependency to 0.13 and removed the `typenum` dependency. [#323](https://github.com/feos-org/feos/pull/323) +- Updated `quantity` dependency to 0.13 and removed the `typenum` dependency. [#328](https://github.com/feos-org/feos/pull/328) ## [Unreleased] ### Added diff --git a/crates/feos-core/src/equation_of_state/mod.rs b/crates/feos-core/src/equation_of_state/mod.rs index dc849068f..14e45c07a 100644 --- a/crates/feos-core/src/equation_of_state/mod.rs +++ b/crates/feos-core/src/equation_of_state/mod.rs @@ -1,7 +1,7 @@ use crate::ReferenceSystem; use crate::state::StateHD; use nalgebra::{ - Const, DVector, DefaultAllocator, Dim, Dyn, OVector, SVector, U1, allocator::Allocator, + Const, DVector, DefaultAllocator, Dim, Dyn, OVector, SVector, allocator::Allocator, }; use num_dual::DualNum; use quantity::{Dimensionless, MolarEnergy, MolarVolume, Temperature}; @@ -114,11 +114,28 @@ impl, D>, D: DualNum + Copy, const N: usize> } /// Ideal gas Helmholtz energy contribution. -pub trait IdealGas { +pub trait IdealGas { /// Implementation of an ideal gas model in terms of the /// logarithm of the cubic thermal de Broglie wavelength /// in units ln(A³) for each component in the system. - fn ln_lambda3 + Copy>(&self, temperature: D2) -> D2; + fn ln_lambda3 + Copy>(&self, temperature: D) -> D; + + /// The name of the ideal gas model. + fn ideal_gas_model(&self) -> &'static str; +} + +/// Ideal gas Helmholtz energy contribution with automatic differentiation with +/// respect to parameters. +pub trait IdealGasAD: Clone { + type Real: IdealGasAD; + type Lifted + Copy>: IdealGasAD; + fn re(&self) -> Self::Real; + fn lift + Copy>(&self) -> Self::Lifted; + + /// Implementation of an ideal gas model in terms of the + /// logarithm of the cubic thermal de Broglie wavelength + /// in units ln(A³) for each component in the system. + fn ln_lambda3(&self, temperature: D) -> D; /// The name of the ideal gas model. fn ideal_gas_model(&self) -> &'static str; @@ -129,45 +146,44 @@ pub trait Total + Copy = f64>: Residual where DefaultAllocator: Allocator, { - type IdealGas: IdealGas; + type RealTotal: Total; + type LiftedTotal + Copy>: Total; + fn re_total(&self) -> Self::RealTotal; + fn lift_total + Copy>(&self) -> Self::LiftedTotal; fn ideal_gas_model(&self) -> &'static str; - fn ideal_gas(&self) -> impl Iterator; - - fn ln_lambda3 + Copy>(&self, temperature: D2) -> OVector { - OVector::from_iterator_generic( - N::from_usize(self.components()), - U1, - self.ideal_gas().map(|i| i.ln_lambda3(temperature)), - ) - } + fn ln_lambda3(&self, temperature: D) -> OVector; - fn ideal_gas_molar_helmholtz_energy + Copy>( + fn ideal_gas_molar_helmholtz_energy( &self, - temperature: D2, - molar_volume: D2, - molefracs: &OVector, - ) -> D2 { + temperature: D, + molar_volume: D, + molefracs: &OVector, + ) -> D { let partial_density = molefracs / molar_volume; - let mut res = D2::from(0.0); - for (i, &r) in self.ideal_gas().zip(partial_density.iter()) { + let mut res = D::from(0.0); + for (&l, &r) in self + .ln_lambda3(temperature) + .iter() + .zip(partial_density.iter()) + { let ln_rho_m1 = if r.re() == 0.0 { - D2::from(0.0) + D::from(0.0) } else { r.ln() - 1.0 }; - res += r * (i.ln_lambda3(temperature) + ln_rho_m1) + res += r * (l + ln_rho_m1) } res * molar_volume * temperature } - fn ideal_gas_helmholtz_energy + Copy>( + fn ideal_gas_helmholtz_energy( &self, - temperature: Temperature, - volume: MolarVolume, - moles: &OVector, - ) -> MolarEnergy { + temperature: Temperature, + volume: MolarVolume, + moles: &OVector, + ) -> MolarEnergy { let total_moles = moles.sum(); let molefracs = moles / total_moles; let molar_volume = volume.into_reduced() / total_moles; @@ -180,32 +196,59 @@ where } impl< - I: IdealGas + Clone + 'static, + I: IdealGas + 'static, C: Deref, R>> + Clone, R: ResidualDyn + 'static, -> Total for C + D: DualNum + Copy, +> Total for C { - type IdealGas = I; + type RealTotal = Self; + type LiftedTotal + Copy> = Self; + fn re_total(&self) -> Self::RealTotal { + self.clone() + } + fn lift_total + Copy>(&self) -> Self::LiftedTotal { + self.clone() + } fn ideal_gas_model(&self) -> &'static str { self.ideal_gas[0].ideal_gas_model() } - fn ideal_gas(&self) -> impl Iterator { - self.ideal_gas.iter() + fn ln_lambda3(&self, temperature: D) -> DVector { + DVector::from_vec( + self.ideal_gas + .iter() + .map(|i| i.ln_lambda3(temperature)) + .collect(), + ) } } -impl + Clone, R: Residual, D>, D: DualNum + Copy, const N: usize> +impl, R: Residual, D>, D: DualNum + Copy, const N: usize> Total, D> for EquationOfState<[I; N], R> { - type IdealGas = I; + type RealTotal = EquationOfState<[I::Real; N], R::Real>; + type LiftedTotal + Copy> = + EquationOfState<[I::Lifted; N], R::Lifted>; + fn re_total(&self) -> Self::RealTotal { + EquationOfState::new( + self.ideal_gas.each_ref().map(|i| i.re()), + self.residual.re(), + ) + } + fn lift_total + Copy>(&self) -> Self::LiftedTotal { + EquationOfState::new( + self.ideal_gas.each_ref().map(|i| i.lift()), + self.residual.lift(), + ) + } fn ideal_gas_model(&self) -> &'static str { self.ideal_gas[0].ideal_gas_model() } - fn ideal_gas(&self) -> impl Iterator { - self.ideal_gas.iter() + fn ln_lambda3(&self, temperature: D) -> SVector { + SVector::from(self.ideal_gas.each_ref().map(|i| i.ln_lambda3(temperature))) } } diff --git a/crates/feos-core/src/lib.rs b/crates/feos-core/src/lib.rs index 48b8306f3..cc280cf0a 100644 --- a/crates/feos-core/src/lib.rs +++ b/crates/feos-core/src/lib.rs @@ -34,8 +34,8 @@ mod phase_equilibria; mod state; pub use ad::{ParametersAD, PropertiesAD}; pub use equation_of_state::{ - EntropyScaling, EquationOfState, IdealGas, Molarweight, NoResidual, Residual, ResidualDyn, - Subset, Total, + EntropyScaling, EquationOfState, IdealGas, IdealGasAD, Molarweight, NoResidual, Residual, + ResidualDyn, Subset, Total, }; pub use errors::{FeosError, FeosResult}; #[cfg(feature = "ndarray")] diff --git a/crates/feos-core/src/phase_equilibria/mod.rs b/crates/feos-core/src/phase_equilibria/mod.rs index 461a1a901..5ed15a767 100644 --- a/crates/feos-core/src/phase_equilibria/mod.rs +++ b/crates/feos-core/src/phase_equilibria/mod.rs @@ -10,7 +10,15 @@ use quantity::{Dimensionless, Energy, Entropy, MolarEnergy, MolarEntropy, Moles} use std::fmt; use std::fmt::Write; +// with empty lines to not mess up the order in the documentation +mod vle_pure; + mod bubble_dew; + +mod tp_flash; + +mod px_flashes; + #[cfg(feature = "ndarray")] mod phase_diagram_binary; #[cfg(feature = "ndarray")] @@ -18,8 +26,7 @@ mod phase_diagram_pure; #[cfg(feature = "ndarray")] mod phase_envelope; mod stability_analysis; -mod tp_flash; -mod vle_pure; + pub use bubble_dew::TemperatureOrPressure; #[cfg(feature = "ndarray")] pub use phase_diagram_binary::PhaseDiagramHetero; @@ -33,22 +40,25 @@ pub use phase_diagram_pure::PhaseDiagram; /// /// ## Contents /// +/// + [Pure component phase equilibria](#pure-component-phase-equilibria) /// + [Bubble and dew point calculations](#bubble-and-dew-point-calculations) -/// + [Heteroazeotropes](#heteroazeotropes) /// + [Flash calculations](#flash-calculations) -/// + [Pure component phase equilibria](#pure-component-phase-equilibria) +/// + [Heteroazeotropes](#heteroazeotropes) /// + [Utility functions](#utility-functions) #[derive(Debug, Clone)] pub struct PhaseEquilibrium + Copy = f64> where DefaultAllocator: Allocator, { - states: [State; P], + pub states: [State; P], pub phase_fractions: [D; P], total_moles: Option>, } -impl fmt::Display for PhaseEquilibrium { +impl, N: Dim, const P: usize> fmt::Display for PhaseEquilibrium +where + DefaultAllocator: Allocator, +{ fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { for (i, s) in self.states.iter().enumerate() { writeln!(f, "phase {i}: {s}")?; @@ -125,12 +135,12 @@ impl, N: Dim, D: DualNum + Copy> PhaseEquilibrium, { - pub(super) fn single_phase(state: State) -> Self { + pub fn single_phase(state: State) -> Self { let total_moles = state.total_moles; Self::with_vapor_phase_fraction(state.clone(), state, D::from(1.0), total_moles) } - pub(super) fn two_phase(vapor: State, liquid: State) -> Self { + pub fn two_phase(vapor: State, liquid: State) -> Self { let (beta, total_moles) = if let (Some(nv), Some(nl)) = (vapor.total_moles, liquid.total_moles) { (nv.convert_into(nl + nv), Some(nl + nv)) @@ -140,7 +150,7 @@ where Self::with_vapor_phase_fraction(vapor, liquid, beta, total_moles) } - pub(super) fn with_vapor_phase_fraction( + pub fn with_vapor_phase_fraction( vapor: State, liquid: State, vapor_phase_fraction: D, @@ -158,11 +168,7 @@ impl, N: Dim, D: DualNum + Copy> PhaseEquilibrium, { - pub(super) fn new( - vapor: State, - liquid1: State, - liquid2: State, - ) -> Self { + pub fn new(vapor: State, liquid1: State, liquid2: State) -> Self { Self { states: [vapor, liquid1, liquid2], phase_fractions: [D::from(1.0), D::from(0.0), D::from(0.0)], @@ -171,7 +177,7 @@ where } } -impl, N: Gradients, const P: usize, D: DualNum + Copy> +impl, N: Gradients, const P: usize, D: DualNum + Copy> PhaseEquilibrium where DefaultAllocator: Allocator, @@ -179,7 +185,13 @@ where pub fn total_moles(&self) -> FeosResult> { self.total_moles.ok_or(FeosError::IntensiveState) } +} +impl, N: Gradients, const P: usize, D: DualNum + Copy> + PhaseEquilibrium +where + DefaultAllocator: Allocator, +{ pub fn molar_enthalpy(&self) -> MolarEnergy { self.states .iter() diff --git a/crates/feos-core/src/phase_equilibria/px_flashes.rs b/crates/feos-core/src/phase_equilibria/px_flashes.rs new file mode 100644 index 000000000..aadeffb76 --- /dev/null +++ b/crates/feos-core/src/phase_equilibria/px_flashes.rs @@ -0,0 +1,465 @@ +#![expect(clippy::toplevel_ref_arg)] +use super::PhaseEquilibrium; +use crate::errors::FeosResult; +use crate::state::State; +use crate::{Composition, FeosError, ReferenceSystem, SolverOptions, Total, Verbosity}; +use nalgebra::allocator::Allocator; +use nalgebra::{DefaultAllocator, Dim, DimAdd, OVector, U1, U2, U3, stack, vector}; +use num_dual::linalg::LU; +use num_dual::{ + Dual, Dual64, DualNum, DualStruct, Gradients, first_derivative, implicit_derivative_sp, partial, +}; +use quantity::{Density, MolarEnergy, MolarEntropy, Pressure, Quantity, SIUnit, Temperature}; + +const MAX_ITER_PX: usize = 20; +const TOL_PX: f64 = 1e-11; + +type PXVars = >::Output; +type TPVars = >::Output; + +impl, N: Gradients + DimAdd + DimAdd, D: DualNum + Copy> + PhaseEquilibrium +where + DefaultAllocator: Allocator + + Allocator + + Allocator> + + Allocator> + + Allocator, PXVars> + + Allocator> + + Allocator> + + Allocator, TPVars>, + PXVars: Gradients, + TPVars: Gradients, +{ + /// Perform a ph-flash calculation. An initial temperature is required + /// and the system needs to be in the two-phase region at that initial + /// temperature. + /// + /// based on Michelsen's work [State function based flash specifications](https://doi.org/10.1016/S0378-3812(99)00092-8) + pub fn ph_flash>( + eos: &E, + pressure: Pressure, + molar_enthalpy: MolarEnergy, + feed: X, + initial_temperature: Temperature, + options: SolverOptions, + ) -> FeosResult { + PhaseEquilibrium::px_flash( + eos, + pressure, + molar_enthalpy, + feed, + initial_temperature, + options, + ) + } + + /// Perform a ps-flash calculation. An initial temperature is required + /// and the system needs to be in the two-phase region at that initial + /// temperature. + /// + /// based on Michelsen's work [State function based flash specifications](https://doi.org/10.1016/S0378-3812(99)00092-8) + pub fn ps_flash>( + eos: &E, + pressure: Pressure, + molar_entropy: MolarEntropy, + feed: X, + initial_temperature: Temperature, + options: SolverOptions, + ) -> FeosResult { + PhaseEquilibrium::px_flash( + eos, + pressure, + molar_entropy, + feed, + initial_temperature, + options, + ) + } + + // Generic implementation of ph and ps flashes. + fn px_flash, U: PXFlash>( + eos: &E, + pressure: Pressure, + specification: Quantity, + feed: X, + initial_temperature: Temperature, + options: SolverOptions, + ) -> FeosResult + where + Quantity: ReferenceSystem, + Quantity: ReferenceSystem, + { + let (max_iter, tol, verbosity) = options.unwrap_or(MAX_ITER_PX, TOL_PX); + let (molefracs, total_moles) = feed.into_molefracs(eos)?; + + // initialize with a tp flash + let eos_f64 = eos.re_total(); + let vle = PhaseEquilibrium::tp_flash( + &eos_f64, + initial_temperature, + pressure.re(), + molefracs.map(|x| x.re()), + None, + Default::default(), + None, + )?; + + // extract specifications + let p = pressure.into_reduced().re(); + let hs = specification.into_reduced().re(); + let z = molefracs.map(|x| x.re()); + let specs = (p, hs, z.clone()); + + // extract variables + let t = initial_temperature.into_reduced(); + let beta = vle.vapor_phase_fraction(); + let rho_v = vle.vapor().density.into_reduced(); + let rho_l = vle.liquid().partial_density().into_reduced(); + let mut vars = stack![rho_l; vector![t, beta, rho_v]]; + let mut old_res = None; + + log_iter!( + verbosity, + " iter | method | temperature | residual | phase I mole fractions | phase II mole fractions " + ); + log_iter!(verbosity, "{:-<102}", ""); + log_iter!( + verbosity, + " {:4} | | {:9.5} | | {:10.8?} | {:10.8?}", + 0, + Temperature::from_reduced(t), + (&rho_l / rho_l.sum() + (&z - &rho_l / rho_l.sum()) / beta).as_slice(), + (&rho_l / rho_l.sum()).as_slice(), + ); + + // iterate + for k in 0..max_iter { + // always try a Newton step first + let (grad, new_vars) = U::newton_step(&eos_f64, &vars, &specs)?; + let new_res = grad.norm(); + let (method, res) = if let Some(r) = old_res + && r < new_res + { + // if the residual is not reduced, reject the step and do a tp-flash instead + vars = U::tp_step(&eos_f64, &vars, &specs)?; + ("Tp-flash", None) + } else { + vars = new_vars; + ("Newton", Some(new_res)) + }; + + if let Verbosity::Iter = verbosity { + let (t, _, _, _, x, y) = unpack_variables(&z, &vars); + log_iter!( + verbosity, + " {:4} | {:^8} | {:9.5} | {} | {:10.8?} | {:10.8?}", + k + 1, + method, + Temperature::from_reduced(t), + res.map_or(String::from(" "), |r| format!("{r:14.8e}")), + y.as_slice(), + x.as_slice(), + ); + } + + if let Some(res) = res + && res < tol + { + log_result!( + verbosity, + "px flash: calculation converged in {} step(s)\n", + k + 1 + ); + + // implicit differentiation + let specs = ( + pressure.into_reduced(), + specification.into_reduced(), + molefracs.clone(), + ); + let vars = implicit_derivative_sp( + |variables, specifications| { + U::state_function(&eos.lift_total(), variables, specifications) + }, + vars, + &specs, + ); + let (t, beta, rho_l, rho_v, x, y) = unpack_variables(&molefracs, &vars); + + // store results in PhaseEquilibrium + let liquid = State::new( + eos, + Temperature::from_reduced(t), + Density::from_reduced(rho_l), + x, + )?; + let vapor = State::new( + eos, + Temperature::from_reduced(t), + Density::from_reduced(rho_v), + y, + )?; + return Ok(PhaseEquilibrium::with_vapor_phase_fraction( + vapor, + liquid, + beta, + total_moles, + )); + } + old_res = res; + } + Err(FeosError::NotConverged("px flash".to_owned())) + } +} + +fn unpack_variables + Copy, N: Dim + DimAdd>( + molefracs: &OVector, + variables: &OVector>, +) -> (D, D, D, D, OVector, OVector) +where + DefaultAllocator: Allocator + Allocator>, +{ + let n = molefracs.len(); + let rho_i_l = variables.rows_generic(0, N::from_usize(n)).clone_owned(); + let [[t, beta, rho_v]] = variables.rows_generic(n, U3).clone_owned().data.0; + let rho_l = rho_i_l.sum(); + let x = rho_i_l / rho_l; + let y = &x + (molefracs - &x) / beta; + (t, beta, rho_l, rho_v, x, y) +} + +fn unpack_tp_variables + Copy, N: Dim + DimAdd>( + molefracs: &OVector, + variables: &OVector>, +) -> (D, D, D, OVector, OVector) +where + DefaultAllocator: Allocator + Allocator>, +{ + let n = molefracs.len(); + let rho_i_l = variables.rows_generic(0, N::from_usize(n)).clone_owned(); + let [[beta, rho_v]] = variables.rows_generic(n, U2).clone_owned().data.0; + let rho_l = rho_i_l.sum(); + let x = rho_i_l / rho_l; + let y = &x + (molefracs - &x) / beta; + (beta, rho_l, rho_v, x, y) +} + +trait PXFlash: Sized + Copy { + // potential function for which the flash solution is a saddle point. + fn state_function, N: Dim + DimAdd, D: DualNum + Copy>( + eos: &E, + variables: OVector>, + args: &(D, D, OVector), + ) -> D + where + DefaultAllocator: Allocator + Allocator>; + + fn evaluate_property, N: Gradients, D: DualNum + Copy>( + vle: &PhaseEquilibrium, + ) -> Quantity + where + DefaultAllocator: Allocator; + + // the potential function for a tp-flash specification (Q = A + V*p_spec) + fn tp_state_function, N: Dim + DimAdd, D: DualNum + Copy>( + eos: &E, + variables: OVector>, + &(t, p, ref z): &(D, D, OVector), + ) -> D + where + DefaultAllocator: Allocator + Allocator>, + { + let (beta, rho_l, rho_v, x, y) = unpack_tp_variables(z, &variables); + let potential = |molefracs, rho: D, t| { + let v = rho.recip(); + let a_res = eos.residual_helmholtz_energy(t, v, &molefracs); + let a_ig = eos.ideal_gas_molar_helmholtz_energy(t, v, &molefracs); + a_res + a_ig + v * p + }; + potential(y, rho_v, t) * beta + potential(x, rho_l, t) * (-beta + 1.0) + } + + // An undamped Newton step for the gradients of the potential function. + // Because the ps and ph flashes are saddle points rather then extrema, + // the value of the potential can not be used as convergence criterion. + #[expect(clippy::type_complexity)] + fn newton_step, N: Dim + DimAdd, D: DualNum + Copy>( + eos: &E, + variables: &OVector>, + specifications: &(D, D, OVector), + ) -> FeosResult<(OVector>, OVector>)> + where + DefaultAllocator: Allocator + Allocator> + Allocator, PXVars>, + PXVars: Gradients, + { + let (_, grad, hess) = PXVars::::hessian( + |variables, specifications| { + Self::state_function(&eos.lift_total(), variables, specifications) + }, + variables, + specifications, + ); + let dx = LU::new(hess)?.solve(&grad); + Ok((grad, variables - &dx)) + } + + // A much slower but more robust step that calculates the implicit + // derivative of the temperature only (which is well behaved + // according to Michelsen) and then calculates all other variables + // from a tp-flash. + fn tp_step, N: Gradients + DimAdd + DimAdd>( + eos: &E, + variables: &OVector>, + &(p, hs_spec, ref z): &(f64, f64, OVector), + ) -> FeosResult>> + where + Quantity: ReferenceSystem, + DefaultAllocator: Allocator + + Allocator + + Allocator> + + Allocator> + + Allocator> + + Allocator, TPVars>, + TPVars: Gradients, + { + let (mut t, beta, rho_l, rho_v, x, y) = unpack_variables(z, variables); + let rho_i_l = rho_l * x; + let (hs, dhs) = first_derivative( + partial( + |t: Dual<_, _>, args: &(_, OVector<_, _>)| { + let &(p, ref z) = args; + let args = (t, p, z.clone_owned()); + + // implicit differentiation of the tp stationarity condition + // to obtain the derivative of the other variables w.r.t. t + let tp_vars = implicit_derivative_sp( + |variables, args| { + Self::tp_state_function(&eos.lift_total().lift_total(), variables, args) + }, + stack![rho_i_l; vector![beta, rho_v]], + &args, + ); + let (beta, rho_l, rho_v, x, y) = unpack_tp_variables(z, &tp_vars); + + // Evaluation of the enthalpy/entropy including the derivatives. + let liquid = State::new( + &eos.lift_total(), + Temperature::from_reduced(t), + Density::from_reduced(rho_l), + x, + )?; + let vapor = State::new( + &eos.lift_total(), + Temperature::from_reduced(t), + Density::from_reduced(rho_v), + y, + )?; + Ok::<_, FeosError>( + Self::evaluate_property(&PhaseEquilibrium::with_vapor_phase_fraction( + vapor, liquid, beta, None, + )) + .into_reduced(), + ) + }, + &(p, z.clone_owned()), + ), + t, + )?; + + // Newton step for the temperature + t -= (hs - hs_spec) / dhs; + + // pack variables into PhaseEquilibrium for initial values + let liquid = State::new_density( + eos, + Temperature::from_reduced(t), + Density::from_reduced(rho_i_l), + )?; + let vapor = State::new( + eos, + Temperature::from_reduced(t), + Density::from_reduced(rho_v), + y, + )?; + let vle = PhaseEquilibrium::with_vapor_phase_fraction(vapor, liquid, beta, None); + + // tp-flash for all other variables + let vle = PhaseEquilibrium::tp_flash( + eos, + Temperature::from_reduced(t), + Pressure::from_reduced(p), + z, + Some(&vle), + Default::default(), + None, + )?; + let beta = vle.vapor_phase_fraction(); + let rho_v = vle.vapor().density.into_reduced(); + let rho_l = vle.liquid().partial_density().into_reduced(); + Ok(stack![rho_l; vector![t, beta, rho_v]]) + } +} + +impl PXFlash for SIUnit<-2, 2, 1, 0, 0, -1, 0> { + // the potential function for a ph-flash specification (Q = (A + V*p_spec - H_spec) / T) + fn state_function, N: Dim + DimAdd, D: DualNum + Copy>( + eos: &E, + variables: OVector>, + &(p, h, ref z): &(D, D, OVector), + ) -> D + where + DefaultAllocator: Allocator + Allocator>, + { + let (t, beta, rho_l, rho_v, x, y) = unpack_variables(z, &variables); + let potential = |molefracs, rho: D, t| { + let v = rho.recip(); + let a_res = eos.residual_helmholtz_energy(t, v, &molefracs); + let a_ig = eos.ideal_gas_molar_helmholtz_energy(t, v, &molefracs); + (a_res + a_ig + v * p - h) / t + }; + potential(y, rho_v, t) * beta + potential(x, rho_l, t) * (-beta + 1.0) + } + + fn evaluate_property, N: Gradients, D: DualNum + Copy>( + vle: &PhaseEquilibrium, + ) -> Quantity + where + DefaultAllocator: Allocator, + { + vle.molar_enthalpy() + } +} + +impl PXFlash for SIUnit<-2, 2, 1, 0, -1, -1, 0> { + // the potential function for a ps-flash specification (Q = A + T*S_spec + V*p_spec) + fn state_function, N: Dim + DimAdd, D: DualNum + Copy>( + eos: &E, + variables: OVector>, + &(p, s, ref z): &(D, D, OVector), + ) -> D + where + DefaultAllocator: Allocator + Allocator>, + { + let (t, beta, rho_l, rho_v, x, y) = unpack_variables(z, &variables); + let potential = |molefracs, rho: D, t| { + let v = rho.recip(); + let a_res = eos.residual_helmholtz_energy(t, v, &molefracs); + let a_ig = eos.ideal_gas_molar_helmholtz_energy(t, v, &molefracs); + // Division by t.re() is done to ensure that the state function has the same + // units (and in conclusion same order of magnitude) as the ph state function. + // This allows using the same toelrances for both methods. + (a_res + a_ig + t * s + v * p) / t.re() + }; + potential(y, rho_v, t) * beta + potential(x, rho_l, t) * (-beta + 1.0) + } + + fn evaluate_property, N: Gradients, D: DualNum + Copy>( + vle: &PhaseEquilibrium, + ) -> Quantity + where + DefaultAllocator: Allocator, + { + vle.molar_entropy() + } +} diff --git a/crates/feos-core/src/state/composition.rs b/crates/feos-core/src/state/composition.rs index f6a190e33..973af8ffb 100644 --- a/crates/feos-core/src/state/composition.rs +++ b/crates/feos-core/src/state/composition.rs @@ -1,4 +1,3 @@ -use super::State; use crate::equation_of_state::Residual; use crate::{FeosError, FeosResult}; use nalgebra::allocator::Allocator; @@ -6,10 +5,27 @@ use nalgebra::{DefaultAllocator, Dim, Dyn, OVector, U1, U2, dvector, vector}; use num_dual::{DualNum, DualStruct}; use quantity::Moles; +/// Trait to generalize over different input types for the composition of +/// a state. +/// +/// The trait is implemented for the following data types: +/// +/// |components|input|total_moles?|comment| +/// |:-:|-|-|-| +/// |1|`()`|-|| +/// |1|`Moles`|✅| +/// |2|`f64`|-| +/// |N|`OVector`|-| +/// |N|`&OVector`|-| +/// |N|`OVector`|-|`Dyn` only| +/// |N|`&OVector`|-|`Dyn` only| +/// |N|`Moles>`|✅| +/// |N|`&Moles>`|✅| pub trait Composition + Copy, N: Dim> where DefaultAllocator: Allocator, { + /// Convert the composition into molefracs and total moles if possible. #[expect(clippy::type_complexity)] fn into_molefracs>( self, @@ -42,19 +58,6 @@ where } } -// copy the composition from a given state -impl + Copy, N: Dim> Composition for &State -where - DefaultAllocator: Allocator, -{ - fn into_molefracs>( - self, - _: &E1, - ) -> FeosResult<(OVector, Option>)> { - Ok(((self.molefracs.clone()), self.total_moles)) - } -} - // a pure component needs no specification impl + Copy> Composition for () { fn into_molefracs>( diff --git a/crates/feos-core/src/state/mod.rs b/crates/feos-core/src/state/mod.rs index 716468aa8..daa7a7ff4 100644 --- a/crates/feos-core/src/state/mod.rs +++ b/crates/feos-core/src/state/mod.rs @@ -306,17 +306,9 @@ where /// Return a new `State` for the combination of inputs. /// - /// The function attempts to create a new state using the given input values. If the state - /// is overdetermined, it will choose a method based on the following hierarchy. - /// 1. Create a state non-iteratively from the set of $T$, $V$, $\rho$, $\rho_i$, $N$, $N_i$ and $x_i$. - /// 2. Use a density iteration for a given pressure. - /// - /// The [StateBuilder] provides a convenient way of calling this function without the need to provide - /// all the optional input values. - /// /// # Errors /// - /// When the state cannot be created using the combination of inputs. + /// When the state cannot be created using the combination of inputs (is over- or underdetermined). pub fn build>( eos: &E, temperature: Temperature, @@ -417,18 +409,9 @@ where { /// Return a new `State` for the combination of inputs. /// - /// The function attempts to create a new state using the given input values. If the state - /// is overdetermined, it will choose a method based on the following hierarchy. - /// 1. Create a state non-iteratively from the set of $T$, $V$, $\rho$, $\rho_i$, $N$, $N_i$ and $x_i$. - /// 2. Use a density iteration for a given pressure. - /// 3. Determine the state using a Newton iteration from (in this order): $(p, h)$, $(p, s)$, $(T, h)$, $(T, s)$, $(V, u)$ - /// - /// The [StateBuilder] provides a convenient way of calling this function without the need to provide - /// all the optional input values. - /// /// # Errors /// - /// When the state cannot be created using the combination of inputs. + /// When the state cannot be created using the combination of inputs (is over- or underdetermined). #[expect(clippy::too_many_arguments)] pub fn build_full + Clone>( eos: &E, @@ -461,28 +444,33 @@ where match state { Some(state) => Ok(state), None => { - // Check if new state can be created using molar_enthalpy and temperature - if let (Some(p), Some(h)) = (pressure, molar_enthalpy) { - return State::new_nph(eos, p, h, composition, density_initialization, ti); - } - if let (Some(p), Some(s)) = (pressure, molar_entropy) { - return State::new_nps(eos, p, s, composition, density_initialization, ti); - } - if let (Some(t), Some(h)) = (temperature, molar_enthalpy) { - return State::new_nth(eos, t, h, composition, density_initialization); - } - if let (Some(t), Some(s)) = (temperature, molar_entropy) { - return State::new_nts(eos, t, s, composition, density_initialization); - } - if let (Some(u), Some(v)) = (molar_internal_energy, volume) { - let (molefracs, total_moles) = composition.into_molefracs(eos)?; - if let Some(n) = total_moles { - return State::new_nvu(eos, v, u, (molefracs, n), ti); + match ( + temperature, + pressure, + volume, + molar_enthalpy, + molar_entropy, + molar_internal_energy, + ) { + (Some(t), None, None, Some(h), None, None) => { + State::new_nth(eos, t, h, composition, density_initialization) + } + (Some(t), None, None, None, Some(s), None) => { + State::new_nts(eos, t, s, composition, density_initialization) + } + (None, Some(p), None, Some(h), None, None) => { + State::new_nph(eos, p, h, composition, density_initialization, ti) + } + (None, Some(p), None, None, Some(s), None) => { + State::new_nps(eos, p, s, composition, density_initialization, ti) + } + (None, None, Some(v), None, None, Some(u)) => { + State::new_nvu(eos, v, u, composition, ti) } + _ => Err(FeosError::UndeterminedState(String::from( + "Missing input parameters.", + ))), } - Err(FeosError::UndeterminedState(String::from( - "Missing input parameters.", - ))) } } } diff --git a/crates/feos-core/src/state/properties.rs b/crates/feos-core/src/state/properties.rs index fa0f53b4a..132bad936 100644 --- a/crates/feos-core/src/state/properties.rs +++ b/crates/feos-core/src/state/properties.rs @@ -20,7 +20,9 @@ where let ideal_gas = || { quantity::ad::gradient_copy( partial2( - |n: Dimensionless<_>, &t, &v| self.eos.ideal_gas_helmholtz_energy(t, v, &n), + |n: Dimensionless<_>, &t, &v| { + self.eos.lift_total().ideal_gas_helmholtz_energy(t, v, &n) + }, &self.temperature, &self.molar_volume, ), @@ -38,7 +40,7 @@ where quantity::ad::partial_hessian_copy( partial( |(n, t): (Dimensionless<_>, _), &v| { - self.eos.ideal_gas_helmholtz_energy(t, v, &n) + self.eos.lift_total().ideal_gas_helmholtz_energy(t, v, &n) }, &self.molar_volume, ), @@ -89,7 +91,7 @@ where let ideal_gas = || { -quantity::ad::first_derivative( partial2( - |t, &v, n| self.eos.ideal_gas_helmholtz_energy(t, v, n), + |t, &v, n| self.eos.lift_total().ideal_gas_helmholtz_energy(t, v, n), &self.molar_volume, &self.molefracs, ), @@ -115,7 +117,7 @@ where let ideal_gas = || { -quantity::ad::second_derivative( partial2( - |t, &v, n| self.eos.ideal_gas_helmholtz_energy(t, v, n), + |t, &v, n| self.eos.lift_total().ideal_gas_helmholtz_energy(t, v, n), &self.molar_volume, &self.molefracs, ), @@ -135,7 +137,7 @@ where let ideal_gas = || { -quantity::ad::third_derivative( partial2( - |t, &v, n| self.eos.ideal_gas_helmholtz_energy(t, v, n), + |t, &v, n| self.eos.lift_total().ideal_gas_helmholtz_energy(t, v, n), &self.molar_volume, &self.molefracs, ), @@ -176,7 +178,7 @@ where let ideal_gas = || { quantity::ad::zeroth_derivative( partial2( - |t, &v, n| self.eos.ideal_gas_helmholtz_energy(t, v, n), + |t, &v, n| self.eos.lift_total().ideal_gas_helmholtz_energy(t, v, n), &self.molar_volume, &self.molefracs, ), @@ -253,7 +255,9 @@ where if let Contributions::IdealGas | Contributions::Total = contributions { res.push(( self.eos.ideal_gas_model(), - self.eos.ideal_gas_molar_helmholtz_energy(t, v, &x), + self.eos + .lift_total() + .ideal_gas_molar_helmholtz_energy(t, v, &x), )); } if let Contributions::Residual | Contributions::Total = contributions { diff --git a/crates/feos-derive/src/ideal_gas.rs b/crates/feos-derive/src/ideal_gas.rs index 06196e82c..f9b5eb671 100644 --- a/crates/feos-derive/src/ideal_gas.rs +++ b/crates/feos-derive/src/ideal_gas.rs @@ -42,7 +42,7 @@ fn impl_ideal_gas( }); quote! { impl IdealGas for IdealGasModel { - fn ln_lambda3 + Copy>(&self, temperature: D) -> D { + fn ln_lambda3 + Copy>(&self, temperature: D) -> D { match self { #(#ln_lambda3,)* } diff --git a/crates/feos-dft/src/profile/properties.rs b/crates/feos-dft/src/profile/properties.rs index 6095fb172..d542b27da 100644 --- a/crates/feos-dft/src/profile/properties.rs +++ b/crates/feos-dft/src/profile/properties.rs @@ -184,7 +184,7 @@ where temperature: Dual64, density: &Array, ) -> Array { - let lambda = self.bulk.eos.ln_lambda3(temperature); + let lambda = self.bulk.eos.lift_total().ln_lambda3(temperature); let mut phi = Array::zeros(density.raw_dim().remove_axis(Axis(0))); for (i, rhoi) in density.outer_iter().enumerate() { phi += &rhoi.mapv(|rhoi| (lambda[i] + rhoi.ln() - 1.0) * rhoi); diff --git a/crates/feos/src/ideal_gas/joback.rs b/crates/feos/src/ideal_gas/joback.rs index 8f2279d5b..e30c630b2 100644 --- a/crates/feos/src/ideal_gas/joback.rs +++ b/crates/feos/src/ideal_gas/joback.rs @@ -1,12 +1,13 @@ //! Implementation of the ideal gas heat capacity (de Broglie wavelength) //! of [Joback and Reid, 1987](https://doi.org/10.1080/00986448708960487). use feos_core::parameter::{FromSegments, Parameters}; -use feos_core::{FeosResult, IdealGas, ReferenceSystem}; +use feos_core::{FeosResult, IdealGas, IdealGasAD, ReferenceSystem}; use nalgebra::DVector; use num_dual::*; use quantity::{MolarEntropy, Temperature}; use serde::{Deserialize, Serialize}; use std::collections::HashMap; +use std::ops::Mul; /// Coefficients used in the Joback model. /// @@ -96,9 +97,9 @@ impl Joback { } } -impl + Copy> IdealGas for Joback { - fn ln_lambda3 + Copy>(&self, temperature: D2) -> D2 { - let [a, b, c, d, e] = self.0.each_ref().map(D2::from_inner); +impl + Copy> Joback { + fn ln_lambda3 + Copy + Mul>(&self, temperature: D2) -> D2 { + let [a, b, c, d, e] = self.0; let t = temperature; let t2 = t * t; let t4 = t2 * t2; @@ -115,6 +116,31 @@ impl + Copy> IdealGas for Joback { + (t / T0).ln() * a; (h - t * s) / (t * RGAS) + f } +} + +impl IdealGas for Joback { + fn ln_lambda3 + Copy>(&self, temperature: D) -> D { + self.ln_lambda3(temperature) + } + + fn ideal_gas_model(&self) -> &'static str { + "Ideal gas (Joback)" + } +} + +impl + Copy> IdealGasAD for Joback { + type Real = Joback; + type Lifted + Copy> = Joback; + fn re(&self) -> Self::Real { + Joback(self.0.each_ref().map(D::re)) + } + fn lift + Copy>(&self) -> Self::Lifted { + Joback(self.0.each_ref().map(D2::from_inner)) + } + + fn ln_lambda3(&self, temperature: D) -> D { + self.ln_lambda3(temperature) + } fn ideal_gas_model(&self) -> &'static str { "Ideal gas (Joback)" diff --git a/crates/feos/src/multiparameter/mod.rs b/crates/feos/src/multiparameter/mod.rs index b4a6c4f9f..fe30093ac 100644 --- a/crates/feos/src/multiparameter/mod.rs +++ b/crates/feos/src/multiparameter/mod.rs @@ -118,10 +118,10 @@ impl Subset for MultiParameter { } impl IdealGas for MultiParameterIdealGas { - fn ln_lambda3 + Copy>(&self, temperature: D2) -> D2 { + fn ln_lambda3 + Copy>(&self, temperature: D) -> D { let tau = temperature.recip() * self.tc; // bit of a hack to convert from phi^0 into ln Lambda^3 - let delta = D2::from(E / (6.02214076e-7 * self.rhoc)); + let delta = D::from(E / (6.02214076e-7 * self.rhoc)); self.terms.iter().map(|r| r.evaluate(delta, tau)).sum() } diff --git a/crates/feos/tests/pcsaft/mod.rs b/crates/feos/tests/pcsaft/mod.rs index 6193aca74..732e6c0c5 100644 --- a/crates/feos/tests/pcsaft/mod.rs +++ b/crates/feos/tests/pcsaft/mod.rs @@ -1,6 +1,7 @@ mod critical_point; mod dft; mod properties; +mod px_flashes; mod stability_analysis; mod state_creation_mixture; mod state_creation_pure; diff --git a/crates/feos/tests/pcsaft/px_flashes.rs b/crates/feos/tests/pcsaft/px_flashes.rs new file mode 100644 index 000000000..4cfb179b6 --- /dev/null +++ b/crates/feos/tests/pcsaft/px_flashes.rs @@ -0,0 +1,138 @@ +use approx::assert_relative_eq; +use feos::ideal_gas::Joback; +use feos::pcsaft::PcSaftBinary; +use feos_core::{ + Contributions, EquationOfState, FeosResult, IdealGasAD, ParametersAD, PhaseEquilibrium, + ReferenceSystem, SolverOptions, Verbosity, +}; +use nalgebra::U1; +use num_dual::{DualStruct, DualVec}; +use quantity::*; + +#[test] +fn test_ph_flash() -> FeosResult<()> { + let params = [ + [1.5, 3.4, 180.0, 2.2, 0.03, 2500., 2.0, 1.0], + [4.5, 3.6, 250.0, 1.2, 0.015, 1500., 1.0, 2.0], + ]; + let kij = 0.15; + let pcsaft = PcSaftBinary::new(params, kij); + let joback = [ + Joback([380., 0.0, 0.0, 0.0, 0.0]), + Joback([210., 0.0, 0.0, 0.0, 0.0]), + ]; + let eos = EquationOfState::new(joback.clone(), pcsaft); + let p = 50.0 * BAR; + let t0 = Some(500.0 * KELVIN); + let x = 0.3; + let dew = PhaseEquilibrium::dew_point(&eos, p, x, t0, None, Default::default())?; + let bubble = PhaseEquilibrium::bubble_point(&eos, p, x, t0, None, Default::default())?; + let h = 0.2 * dew.molar_enthalpy() + 0.8 * bubble.molar_enthalpy(); + let t0 = 0.8 * dew.vapor().temperature + 0.2 * bubble.vapor().temperature; + let options = SolverOptions { + verbosity: Verbosity::Iter, + ..Default::default() + }; + let vle = PhaseEquilibrium::ph_flash(&eos, p, h, x, t0, options)?; + println!("{vle}"); + println!("{h}\n{}", vle.molar_enthalpy()); + assert_relative_eq!(h, vle.molar_enthalpy(), max_relative = 1e-10); + + let pcsaft_ad = pcsaft.named_derivatives(["k_ij"]); + let joback_ad = joback.each_ref().map(|j| j.lift()); + let eos_ad = EquationOfState::new(joback_ad, pcsaft_ad); + let vle_ad = PhaseEquilibrium::ph_flash( + &eos_ad, + Pressure::from_inner(&p), + MolarEnergy::from_inner(&h), + DualVec::from_inner(&x), + t0, + Default::default(), + )?; + let [[dt]] = vle_ad + .vapor() + .temperature + .into_reduced() + .eps + .unwrap_generic(U1, U1) + .data + .0; + println!("{dt}"); + + let dkij = 1e-7; + let pcsaft_h = PcSaftBinary::new(params, kij + dkij); + let eos_h = EquationOfState::new(joback.clone(), pcsaft_h); + let vle_h = PhaseEquilibrium::ph_flash(&eos_h, p, h, x, t0, Default::default())?; + let dt_h = (vle_h.vapor().temperature - vle.vapor().temperature).into_reduced() / dkij; + println!("{dt_h}"); + assert_relative_eq!(dt, dt_h, max_relative = 1e-4); + + Ok(()) +} + +#[test] +fn test_ps_flash() -> FeosResult<()> { + let params = [ + [1.5, 3.4, 180.0, 2.2, 0.03, 2500., 2.0, 1.0], + [4.5, 3.6, 250.0, 1.2, 0.015, 1500., 1.0, 2.0], + ]; + let kij = 0.15; + let pcsaft = PcSaftBinary::new(params, kij); + let joback = [ + Joback([380., 0.0, 0.0, 0.0, 0.0]), + Joback([210., 0.0, 0.0, 0.0, 0.0]), + ]; + let eos = EquationOfState::new(joback.clone(), pcsaft); + let p = 50.0 * BAR; + println!( + "{}", + PhaseEquilibrium::bubble_point(&eos, 500.0 * KELVIN, 0.5, None, None, Default::default())? + .vapor() + .pressure(Contributions::Total) + ); + let t0 = Some(500.0 * KELVIN); + let x = 0.3; + let dew = PhaseEquilibrium::dew_point(&eos, p, x, t0, None, Default::default())?; + let bubble = PhaseEquilibrium::bubble_point(&eos, p, x, t0, None, Default::default())?; + let s = 0.2 * dew.molar_entropy() + 0.8 * bubble.molar_entropy(); + let t0 = 0.8 * dew.vapor().temperature + 0.2 * bubble.vapor().temperature; + let options = SolverOptions { + verbosity: Verbosity::Iter, + ..Default::default() + }; + let vle = PhaseEquilibrium::ps_flash(&eos, p, s, x, t0, options)?; + println!("{vle}"); + println!("{s}\n{}", vle.molar_entropy()); + assert_relative_eq!(s, vle.molar_entropy(), max_relative = 1e-10); + + let pcsaft_ad = pcsaft.named_derivatives(["k_ij"]); + let joback_ad = joback.each_ref().map(|j| j.lift()); + let eos_ad = EquationOfState::new(joback_ad, pcsaft_ad); + let vle_ad = PhaseEquilibrium::ps_flash( + &eos_ad, + Pressure::from_inner(&p), + MolarEntropy::from_inner(&s), + DualVec::from_inner(&x), + t0, + Default::default(), + )?; + let [[dt]] = vle_ad + .vapor() + .temperature + .into_reduced() + .eps + .unwrap_generic(U1, U1) + .data + .0; + println!("{dt}"); + + let dkij = 1e-7; + let pcsaft_h = PcSaftBinary::new(params, kij + dkij); + let eos_h = EquationOfState::new(joback.clone(), pcsaft_h); + let vle_h = PhaseEquilibrium::ps_flash(&eos_h, p, s, x, t0, Default::default())?; + let dt_h = (vle_h.vapor().temperature - vle.vapor().temperature).into_reduced() / dkij; + println!("{dt_h}"); + assert_relative_eq!(dt, dt_h, max_relative = 1e-4); + + Ok(()) +} diff --git a/crates/feos/tests/pcsaft/tp_flash.rs b/crates/feos/tests/pcsaft/tp_flash.rs index b31a14998..0f28e0c2e 100644 --- a/crates/feos/tests/pcsaft/tp_flash.rs +++ b/crates/feos/tests/pcsaft/tp_flash.rs @@ -32,15 +32,7 @@ fn test_tp_flash() -> FeosResult<()> { println!("{p_propane} {p_butane} {x1} {y1} {z1}"); let mix = PcSaft::new(read_params(vec!["propane", "butane"])?); let options = SolverOptions::new().max_iter(100).tol(1e-12); - let vle = PhaseEquilibrium::tp_flash( - &&mix, - t, - p, - &(dvector![z1, 1.0 - z1] * MOL), - None, - options, - None, - )?; + let vle = PhaseEquilibrium::tp_flash(&&mix, t, p, z1, None, options, None)?; println!( "x1: {}, y1: {}", vle.liquid().molefracs[0], @@ -85,7 +77,7 @@ fn test_tp_flash_zero_component() -> FeosResult<()> { &&eos_full, 300.0 * KELVIN, 1.2 * BAR, - &(dvector![0.0, 0.5, 0.5] * MOL), + dvector![0.0, 0.5, 0.5], None, options, None, @@ -94,7 +86,7 @@ fn test_tp_flash_zero_component() -> FeosResult<()> { &&eos_binary, 300.0 * KELVIN, 1.2 * BAR, - &(dvector![0.5, 0.5] * MOL), + dvector![0.5, 0.5], None, options, None, diff --git a/docs/recipes/index.md b/docs/recipes/index.md index 5a7a95086..65b3a340a 100644 --- a/docs/recipes/index.md +++ b/docs/recipes/index.md @@ -12,6 +12,7 @@ If you are looking for tutorials with explanations, see the [tutorials](/tutoria recipes_critical_point_pure recipes_p_sat_t_boil recipes_phase_equilibrium_pure + recipes_phase_equilibrium_flash recipes_phase_diagram_pure recipes_automatic_differentiation ``` diff --git a/docs/recipes/recipes_phase_equilibrium_flash.ipynb b/docs/recipes/recipes_phase_equilibrium_flash.ipynb new file mode 100644 index 000000000..60158b47f --- /dev/null +++ b/docs/recipes/recipes_phase_equilibrium_flash.ipynb @@ -0,0 +1,149 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "8bec74cc", + "metadata": {}, + "outputs": [], + "source": [ + "import si_units as si\n", + "import feos\n", + "\n", + "parameters = feos.Parameters.from_json(\n", + " substances=['methanol', '1-propanol'], \n", + " pure_path='../../parameters/pcsaft/gross2002.json'\n", + ")\n", + "ideal_gas_parameters = feos.Parameters.from_json(\n", + " substances=['methanol', '1-propanol'], \n", + " pure_path='../../parameters/ideal_gas/poling2000.json'\n", + ")\n", + "eos = feos.EquationOfState.pcsaft(parameters).dippr(ideal_gas_parameters)" + ] + }, + { + "cell_type": "markdown", + "id": "0ace1cfd", + "metadata": {}, + "source": [ + "## Tp-flash" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e11fa945", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " dew point pressure: 0.6356 bar\n", + "bubble point pressure: 0.9245 bar\n" + ] + }, + { + "data": { + "text/markdown": [ + "||temperature|density|molefracs|\n", + "|-|-|-|-|\n", + "|phase 1|350.00000 K|28.75751 mol/m³|[0.59847, 0.40153]|\n", + "|phase 2|350.00000 K|14.29164 kmol/m³|[0.28954, 0.71046]|\n" + ], + "text/plain": [ + "phase 0: T = 350.00000 K, ρ = 28.75751 mol/m³, x = [0.59847, 0.40153]\n", + "phase 1: T = 350.00000 K, ρ = 14.29164 kmol/m³, x = [0.28954, 0.71046]" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "x1 = 0.4\n", + "temperature = 350*si.KELVIN\n", + "\n", + "p_bubble = feos.PhaseEquilibrium.bubble_point(eos, temperature, x1).liquid.pressure()\n", + "p_dew = feos.PhaseEquilibrium.dew_point(eos, temperature, x1).vapor.pressure()\n", + "print(f\"bubble point pressure: {p_bubble/si.BAR:.4} bar\")\n", + "print(f\" dew point pressure: {p_dew/si.BAR:.4} bar\")\n", + "\n", + "feos.PhaseEquilibrium.tp_flash(eos, temperature, 0.8*si.BAR, x1)" + ] + }, + { + "cell_type": "markdown", + "id": "73c81e37", + "metadata": {}, + "source": [ + "## ph-flash" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b3acf88f", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "bubble point\ttemperature: 352.06 K\tenthalpy: -36.5900 kJ/mol\n", + " dew point\ttemperature: 361.09 K\tenthalpy: 3.9800 kJ/mol\n" + ] + }, + { + "data": { + "text/markdown": [ + "||temperature|density|molefracs|\n", + "|-|-|-|-|\n", + "|phase 1|360.47941 K|34.58586 mol/m³|[0.42362, 0.57638]|\n", + "|phase 2|360.47941 K|13.28059 kmol/m³|[0.17482, 0.82518]|\n" + ], + "text/plain": [ + "phase 0: T = 360.47941 K, ρ = 34.58586 mol/m³, x = [0.42362, 0.57638]\n", + "phase 1: T = 360.47941 K, ρ = 13.28059 kmol/m³, x = [0.17482, 0.82518]" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "x1 = 0.4\n", + "pressure = si.BAR\n", + "bubble = feos.PhaseEquilibrium.bubble_point(eos, pressure, x1, 300*si.KELVIN).liquid\n", + "dew = feos.PhaseEquilibrium.dew_point(eos, pressure, x1, 300*si.KELVIN).vapor\n", + "print(f\"bubble point\\ttemperature: {bubble.temperature/si.KELVIN:.2f} K\\tenthalpy: {bubble.molar_enthalpy()/(si.KILO*si.JOULE/si.MOL):8.4f} kJ/mol\")\n", + "print(f\" dew point\\ttemperature: {dew.temperature/si.KELVIN:.2f} K\\tenthalpy: {dew.molar_enthalpy()/(si.KILO*si.JOULE/si.MOL):8.4f} kJ/mol\")\n", + "\n", + "feos.PhaseEquilibrium.ph_flash(eos, pressure, 0*si.JOULE/si.MOL, x1, 356*si.KELVIN)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "feos_devel", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.14.0" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/py-feos/src/phase_equilibria.rs b/py-feos/src/phase_equilibria.rs index c6c857506..feaa38930 100644 --- a/py-feos/src/phase_equilibria.rs +++ b/py-feos/src/phase_equilibria.rs @@ -159,6 +159,128 @@ impl PyPhaseEquilibrium { )) } + /// Create a liquid and vapor state in equilibrium + /// for given pressure, enthalpy and feed composition. + /// + /// Can also be used to calculate liquid liquid phase separation. + /// + /// Parameters + /// ---------- + /// eos : EquationOfState + /// The equation of state. + /// pressure : SINumber + /// The system pressure. + /// molar_enthalpy : SINumber + /// The molar enthalpy of the system. + /// feed : float | SINumber | numpy.ndarray[float] | SIArray1 | list[float] + /// Feed composition. + /// initial_temperature : SINumber + /// The system temperature. + /// max_iter : int, optional + /// The maximum number of iterations. + /// tol: float, optional + /// The solution tolerance. + /// verbosity : Verbosity, optional + /// The verbosity. + /// + /// Returns + /// ------- + /// PhaseEquilibrium + /// + /// Raises + /// ------ + /// RuntimeError + /// When pressure iteration fails or no phase equilibrium is found. + #[staticmethod] + #[pyo3( + text_signature = "(eos, pressure, molar_enthalpy, feed, initial_temperature, max_iter=None, tol=None, verbosity=None)" + )] + #[pyo3(signature = (eos, pressure, molar_enthalpy, feed, initial_temperature, max_iter=None, tol=None, verbosity=None))] + #[expect(clippy::too_many_arguments)] + pub(crate) fn ph_flash( + eos: &PyEquationOfState, + pressure: Pressure, + molar_enthalpy: MolarEnergy, + feed: &Bound<'_, PyAny>, + initial_temperature: Temperature, + max_iter: Option, + tol: Option, + verbosity: Option, + ) -> PyResult { + Ok(Self( + PhaseEquilibrium::ph_flash( + &eos.0, + pressure, + molar_enthalpy, + Compositions::try_from(Some(feed))?, + initial_temperature, + (max_iter, tol, verbosity.map(|v| v.into())).into(), + ) + .map_err(PyFeosError::from)?, + )) + } + + /// Create a liquid and vapor state in equilibrium + /// for given pressure, entropy and feed composition. + /// + /// Can also be used to calculate liquid liquid phase separation. + /// + /// Parameters + /// ---------- + /// eos : EquationOfState + /// The equation of state. + /// pressure : SINumber + /// The system pressure. + /// molar_entropy : SINumber + /// The molar entropy of the system. + /// feed : float | SINumber | numpy.ndarray[float] | SIArray1 | list[float] + /// Feed composition. + /// initial_temperature : SINumber + /// The system temperature. + /// max_iter : int, optional + /// The maximum number of iterations. + /// tol: float, optional + /// The solution tolerance. + /// verbosity : Verbosity, optional + /// The verbosity. + /// + /// Returns + /// ------- + /// PhaseEquilibrium + /// + /// Raises + /// ------ + /// RuntimeError + /// When pressure iteration fails or no phase equilibrium is found. + #[staticmethod] + #[pyo3( + text_signature = "(eos, pressure, molar_entropy, feed, initial_temperature, max_iter=None, tol=None, verbosity=None)" + )] + #[pyo3(signature = (eos, pressure, molar_entropy, feed, initial_temperature, max_iter=None, tol=None, verbosity=None))] + #[expect(clippy::too_many_arguments)] + pub(crate) fn ps_flash( + eos: &PyEquationOfState, + pressure: Pressure, + molar_entropy: MolarEntropy, + feed: &Bound<'_, PyAny>, + initial_temperature: Temperature, + max_iter: Option, + tol: Option, + verbosity: Option, + ) -> PyResult { + Ok(Self( + PhaseEquilibrium::ps_flash( + &eos.0, + pressure, + molar_entropy, + Compositions::try_from(Some(feed))?, + initial_temperature, + (max_iter, tol, verbosity.map(|v| v.into())).into(), + ) + .map_err(PyFeosError::from)?, + )) + } + /// Compute a phase equilibrium for given temperature /// or pressure and liquid mole fractions. /// @@ -345,6 +467,36 @@ impl PyPhaseEquilibrium { PyState(self.0.liquid().clone()) } + #[getter] + fn get_vapor_phase_fraction(&self) -> f64 { + self.0.vapor_phase_fraction() + } + + #[getter] + fn get_total_moles(&self) -> PyResult { + Ok(self.0.total_moles().map_err(PyFeosError::from)?) + } + + #[getter] + fn get_molar_enthalpy(&self) -> MolarEnergy { + self.0.molar_enthalpy() + } + + #[getter] + fn get_enthalpy(&self) -> PyResult { + Ok(self.0.enthalpy().map_err(PyFeosError::from)?) + } + + #[getter] + fn get_molar_entropy(&self) -> MolarEntropy { + self.0.molar_entropy() + } + + #[getter] + fn get_entropy(&self) -> PyResult { + Ok(self.0.entropy().map_err(PyFeosError::from)?) + } + /// Calculate the pure component vapor-liquid equilibria for all /// components in the system. ///