Implement ph and ps flashes for binary mixtures

feos-org · prehner · Mar 31, 2026 · Mar 15, 2026 · Mar 31, 2026 · 268136fa0235a79e3b3a726b405ead28bb2fd1b8
commit 268136fa0235a79e3b3a726b405ead28bb2fd1b8
diff --git 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
@@ -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<I: Clone, R: Residual<Const<N>, D>, D: DualNum<f64> + Copy, const N: usize>
 }
 
 /// Ideal gas Helmholtz energy contribution.
-pub trait IdealGas<D = f64> {
+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<D2: DualNum<f64, Inner = D> + Copy>(&self, temperature: D2) -> D2;
+    fn ln_lambda3<D: DualNum<f64> + 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<D = f64>: Clone {
+    type Real: IdealGasAD;
+    type Lifted<D2: DualNum<f64, Inner = D> + Copy>: IdealGasAD<D2>;
+    fn re(&self) -> Self::Real;
+    fn lift<D2: DualNum<f64, Inner = D> + Copy>(&self) -> Self::Lifted<D2>;
+
+    /// 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<N: Dim = Dyn, D: DualNum<f64> + Copy = f64>: Residual<N, D>
 where
     DefaultAllocator: Allocator<N>,
 {
-    type IdealGas: IdealGas<D>;
+    type RealTotal: Total<N, f64>;
+    type LiftedTotal<D2: DualNum<f64, Inner = D> + Copy>: Total<N, D2>;
+    fn re_total(&self) -> Self::RealTotal;
+    fn lift_total<D2: DualNum<f64, Inner = D> + Copy>(&self) -> Self::LiftedTotal<D2>;
 
     fn ideal_gas_model(&self) -> &'static str;
 
-    fn ideal_gas(&self) -> impl Iterator<Item = &Self::IdealGas>;
-
-    fn ln_lambda3<D2: DualNum<f64, Inner = D> + Copy>(&self, temperature: D2) -> OVector<D2, N> {
-        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<D, N>;
 
-    fn ideal_gas_molar_helmholtz_energy<D2: DualNum<f64, Inner = D> + Copy>(
+    fn ideal_gas_molar_helmholtz_energy(
         &self,
-        temperature: D2,
-        molar_volume: D2,
-        molefracs: &OVector<D2, N>,
-    ) -> D2 {
+        temperature: D,
+        molar_volume: D,
+        molefracs: &OVector<D, N>,
+    ) -> 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<D2: DualNum<f64, Inner = D> + Copy>(
+    fn ideal_gas_helmholtz_energy(
         &self,
-        temperature: Temperature<D2>,
-        volume: MolarVolume<D2>,
-        moles: &OVector<D2, N>,
-    ) -> MolarEnergy<D2> {
+        temperature: Temperature<D>,
+        volume: MolarVolume<D>,
+        moles: &OVector<D, N>,
+    ) -> MolarEnergy<D> {
         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<Target = EquationOfState<Vec<I>, R>> + Clone,
     R: ResidualDyn + 'static,
-> Total<Dyn, f64> for C
+    D: DualNum<f64> + Copy,
+> Total<Dyn, D> for C
 {
-    type IdealGas = I;
+    type RealTotal = Self;
+    type LiftedTotal<D2: DualNum<f64, Inner = D> + Copy> = Self;
+    fn re_total(&self) -> Self::RealTotal {
+        self.clone()
+    }
+    fn lift_total<D2: DualNum<f64, Inner = D> + Copy>(&self) -> Self::LiftedTotal<D2> {
+        self.clone()
+    }
 
     fn ideal_gas_model(&self) -> &'static str {
         self.ideal_gas[0].ideal_gas_model()
     }
 
-    fn ideal_gas(&self) -> impl Iterator<Item = &Self::IdealGas> {
-        self.ideal_gas.iter()
+    fn ln_lambda3(&self, temperature: D) -> DVector<D> {
+        DVector::from_vec(
+            self.ideal_gas
+                .iter()
+                .map(|i| i.ln_lambda3(temperature))
+                .collect(),
+        )
     }
 }
 
-impl<I: IdealGas<D> + Clone, R: Residual<Const<N>, D>, D: DualNum<f64> + Copy, const N: usize>
+impl<I: IdealGasAD<D>, R: Residual<Const<N>, D>, D: DualNum<f64> + Copy, const N: usize>
     Total<Const<N>, D> for EquationOfState<[I; N], R>
 {
-    type IdealGas = I;
+    type RealTotal = EquationOfState<[I::Real; N], R::Real>;
+    type LiftedTotal<D2: DualNum<f64, Inner = D> + Copy> =
+        EquationOfState<[I::Lifted<D2>; N], R::Lifted<D2>>;
+    fn re_total(&self) -> Self::RealTotal {
+        EquationOfState::new(
+            self.ideal_gas.each_ref().map(|i| i.re()),
+            self.residual.re(),
+        )
+    }
+    fn lift_total<D2: DualNum<f64, Inner = D> + Copy>(&self) -> Self::LiftedTotal<D2> {
+        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<Item = &Self::IdealGas> {
-        self.ideal_gas.iter()
+    fn ln_lambda3(&self, temperature: D) -> SVector<D, N> {
+        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
@@ -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
@@ -10,16 +10,23 @@ 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")]
 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<E, const P: usize, N: Dim = Dyn, D: DualNum<f64> + Copy = f64>
 where
     DefaultAllocator: Allocator<N>,
 {
-    states: [State<E, N, D>; P],
+    pub states: [State<E, N, D>; P],
     pub phase_fractions: [D; P],
     total_moles: Option<Moles<D>>,
 }
 
-impl<E: Residual, const P: usize> fmt::Display for PhaseEquilibrium<E, P> {
+impl<E: Residual<N>, N: Dim, const P: usize> fmt::Display for PhaseEquilibrium<E, P, N>
+where
+    DefaultAllocator: Allocator<N>,
+{
     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<E: Residual<N, D>, N: Dim, D: DualNum<f64> + Copy> PhaseEquilibrium<E, 2, N
 where
     DefaultAllocator: Allocator<N>,
 {
-    pub(super) fn single_phase(state: State<E, N, D>) -> Self {
+    pub fn single_phase(state: State<E, N, D>) -> 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<E, N, D>, liquid: State<E, N, D>) -> Self {
+    pub fn two_phase(vapor: State<E, N, D>, liquid: State<E, N, D>) -> 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<E, N, D>,
         liquid: State<E, N, D>,
         vapor_phase_fraction: D,
@@ -158,11 +168,7 @@ impl<E: Residual<N, D>, N: Dim, D: DualNum<f64> + Copy> PhaseEquilibrium<E, 3, N
 where
     DefaultAllocator: Allocator<N>,
 {
-    pub(super) fn new(
-        vapor: State<E, N, D>,
-        liquid1: State<E, N, D>,
-        liquid2: State<E, N, D>,
-    ) -> Self {
+    pub fn new(vapor: State<E, N, D>, liquid1: State<E, N, D>, liquid2: State<E, N, D>) -> Self {
         Self {
             states: [vapor, liquid1, liquid2],
             phase_fractions: [D::from(1.0), D::from(0.0), D::from(0.0)],
@@ -171,15 +177,21 @@ where
     }
 }
 
-impl<E: Total<N, D>, N: Gradients, const P: usize, D: DualNum<f64> + Copy>
+impl<E: Residual<N, D>, N: Gradients, const P: usize, D: DualNum<f64> + Copy>
     PhaseEquilibrium<E, P, N, D>
 where
     DefaultAllocator: Allocator<N>,
 {
     pub fn total_moles(&self) -> FeosResult<Moles<D>> {
         self.total_moles.ok_or(FeosError::IntensiveState)
     }
+}
 
+impl<E: Total<N, D>, N: Gradients, const P: usize, D: DualNum<f64> + Copy>
+    PhaseEquilibrium<E, P, N, D>
+where
+    DefaultAllocator: Allocator<N>,
+{
     pub fn molar_enthalpy(&self) -> MolarEnergy<D> {
         self.states
             .iter()