feos_core/phase_equilibria/
mod.rs

1use crate::equation_of_state::Residual;
2use crate::errors::{EosError, EosResult};
3use crate::state::{DensityInitialization, State};
4use crate::{Contributions, ReferenceSystem};
5use ndarray::Array1;
6use quantity::{Dimensionless, Energy, Moles, Pressure, Temperature, RGAS};
7use std::fmt;
8use std::fmt::Write;
9use std::sync::Arc;
10
11mod bubble_dew;
12mod phase_diagram_binary;
13mod phase_diagram_pure;
14mod phase_envelope;
15mod stability_analysis;
16mod tp_flash;
17mod vle_pure;
18pub use bubble_dew::TemperatureOrPressure;
19pub use phase_diagram_binary::PhaseDiagramHetero;
20pub use phase_diagram_pure::PhaseDiagram;
21
22/// A thermodynamic equilibrium state.
23///
24/// The struct is parametrized over the number of phases with most features
25/// being implemented for the two phase vapor/liquid or liquid/liquid case.
26///
27/// ## Contents
28///
29/// + [Bubble and dew point calculations](#bubble-and-dew-point-calculations)
30/// + [Heteroazeotropes](#heteroazeotropes)
31/// + [Flash calculations](#flash-calculations)
32/// + [Pure component phase equilibria](#pure-component-phase-equilibria)
33/// + [Utility functions](#utility-functions)
34#[derive(Debug)]
35pub struct PhaseEquilibrium<E, const N: usize>([State<E>; N]);
36
37impl<E, const N: usize> Clone for PhaseEquilibrium<E, N> {
38    fn clone(&self) -> Self {
39        Self(self.0.clone())
40    }
41}
42
43impl<E: Residual, const N: usize> fmt::Display for PhaseEquilibrium<E, N> {
44    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
45        for (i, s) in self.0.iter().enumerate() {
46            writeln!(f, "phase {}: {}", i, s)?;
47        }
48        Ok(())
49    }
50}
51
52impl<E: Residual, const N: usize> PhaseEquilibrium<E, N> {
53    pub fn _repr_markdown_(&self) -> String {
54        if self.0[0].eos.components() == 1 {
55            let mut res = "||temperature|density|\n|-|-|-|\n".to_string();
56            for (i, s) in self.0.iter().enumerate() {
57                writeln!(
58                    res,
59                    "|phase {}|{:.5}|{:.5}|",
60                    i + 1,
61                    s.temperature,
62                    s.density
63                )
64                .unwrap();
65            }
66            res
67        } else {
68            let mut res = "||temperature|density|molefracs|\n|-|-|-|-|\n".to_string();
69            for (i, s) in self.0.iter().enumerate() {
70                writeln!(
71                    res,
72                    "|phase {}|{:.5}|{:.5}|{:.5}|",
73                    i + 1,
74                    s.temperature,
75                    s.density,
76                    s.molefracs
77                )
78                .unwrap();
79            }
80            res
81        }
82    }
83}
84
85impl<E> PhaseEquilibrium<E, 2> {
86    pub fn vapor(&self) -> &State<E> {
87        &self.0[0]
88    }
89
90    pub fn liquid(&self) -> &State<E> {
91        &self.0[1]
92    }
93}
94
95impl<E> PhaseEquilibrium<E, 3> {
96    pub fn vapor(&self) -> &State<E> {
97        &self.0[0]
98    }
99
100    pub fn liquid1(&self) -> &State<E> {
101        &self.0[1]
102    }
103
104    pub fn liquid2(&self) -> &State<E> {
105        &self.0[2]
106    }
107}
108
109impl<E: Residual> PhaseEquilibrium<E, 2> {
110    pub(super) fn from_states(state1: State<E>, state2: State<E>) -> Self {
111        let (vapor, liquid) = if state1.density < state2.density {
112            (state1, state2)
113        } else {
114            (state2, state1)
115        };
116        Self([vapor, liquid])
117    }
118
119    /// Creates a new PhaseEquilibrium that contains two states at the
120    /// specified temperature, pressure and moles.
121    ///
122    /// The constructor can be used in custom phase equilibrium solvers or,
123    /// e.g., to generate initial guesses for an actual VLE solver.
124    /// In general, the two states generated are NOT in an equilibrium.
125    pub fn new_npt(
126        eos: &Arc<E>,
127        temperature: Temperature,
128        pressure: Pressure,
129        vapor_moles: &Moles<Array1<f64>>,
130        liquid_moles: &Moles<Array1<f64>>,
131    ) -> EosResult<Self> {
132        let liquid = State::new_npt(
133            eos,
134            temperature,
135            pressure,
136            liquid_moles,
137            DensityInitialization::Liquid,
138        )?;
139        let vapor = State::new_npt(
140            eos,
141            temperature,
142            pressure,
143            vapor_moles,
144            DensityInitialization::Vapor,
145        )?;
146        Ok(Self([vapor, liquid]))
147    }
148
149    pub(super) fn vapor_phase_fraction(&self) -> f64 {
150        (self.vapor().total_moles / (self.vapor().total_moles + self.liquid().total_moles))
151            .into_value()
152    }
153}
154
155impl<E: Residual, const N: usize> PhaseEquilibrium<E, N> {
156    pub(super) fn update_pressure(
157        mut self,
158        temperature: Temperature,
159        pressure: Pressure,
160    ) -> EosResult<Self> {
161        for s in self.0.iter_mut() {
162            *s = State::new_npt(
163                &s.eos,
164                temperature,
165                pressure,
166                &s.moles,
167                DensityInitialization::InitialDensity(s.density),
168            )?;
169        }
170        Ok(self)
171    }
172
173    pub(super) fn update_moles(
174        &mut self,
175        pressure: Pressure,
176        moles: [&Moles<Array1<f64>>; N],
177    ) -> EosResult<()> {
178        for (i, s) in self.0.iter_mut().enumerate() {
179            *s = State::new_npt(
180                &s.eos,
181                s.temperature,
182                pressure,
183                moles[i],
184                DensityInitialization::InitialDensity(s.density),
185            )?;
186        }
187        Ok(())
188    }
189
190    // Total Gibbs energy excluding the constant contribution RT sum_i N_i ln(\Lambda_i^3)
191    pub(super) fn total_gibbs_energy(&self) -> Energy {
192        self.0.iter().fold(Energy::from_reduced(0.0), |acc, s| {
193            let ln_rho = s.partial_density.to_reduced().mapv(f64::ln);
194            acc + s.residual_helmholtz_energy()
195                + s.pressure(Contributions::Total) * s.volume
196                + RGAS * s.temperature * (s.moles.clone() * Dimensionless::new(ln_rho - 1.0)).sum()
197        })
198    }
199}
200
201const TRIVIAL_REL_DEVIATION: f64 = 1e-5;
202
203/// # Utility functions
204impl<E: Residual> PhaseEquilibrium<E, 2> {
205    pub(super) fn check_trivial_solution(self) -> EosResult<Self> {
206        if Self::is_trivial_solution(self.vapor(), self.liquid()) {
207            Err(EosError::TrivialSolution)
208        } else {
209            Ok(self)
210        }
211    }
212
213    /// Check if the two states form a trivial solution
214    pub fn is_trivial_solution(state1: &State<E>, state2: &State<E>) -> bool {
215        let rho1 = state1.partial_density.to_reduced();
216        let rho2 = state2.partial_density.to_reduced();
217
218        rho1.iter()
219            .zip(rho2.iter())
220            .fold(0.0, |acc, (&rho1, &rho2)| {
221                (rho2 / rho1 - 1.0).abs().max(acc)
222            })
223            < TRIVIAL_REL_DEVIATION
224    }
225}
feos_core/phase_equilibria/mod.rs

feos_core/phase_equilibria/
mod.rs
