feos_core/phase_equilibria/
mod.rs

1use crate::FeosError;
2use crate::equation_of_state::Residual;
3use crate::errors::FeosResult;
4use crate::state::State;
5use crate::{Contributions::Total as Tot, ReferenceSystem, Total};
6use nalgebra::allocator::Allocator;
7use nalgebra::{DefaultAllocator, Dim, Dyn};
8use num_dual::{DualNum, Gradients};
9use quantity::{Dimensionless, Energy, Entropy, MolarEnergy, MolarEntropy, Moles};
10use std::fmt;
11use std::fmt::Write;
12
13// with empty lines to not mess up the order in the documentation
14mod vle_pure;
15
16mod bubble_dew;
17
18mod tp_flash;
19
20mod px_flashes;
21
22#[cfg(feature = "ndarray")]
23mod phase_diagram_binary;
24#[cfg(feature = "ndarray")]
25mod phase_diagram_pure;
26#[cfg(feature = "ndarray")]
27mod phase_envelope;
28mod stability_analysis;
29
30pub use bubble_dew::TemperatureOrPressure;
31#[cfg(feature = "ndarray")]
32pub use phase_diagram_binary::PhaseDiagramHetero;
33#[cfg(feature = "ndarray")]
34pub use phase_diagram_pure::PhaseDiagram;
35
36/// A thermodynamic equilibrium state.
37///
38/// The struct is parametrized over the number of phases with most features
39/// being implemented for the two phase vapor/liquid or liquid/liquid case.
40///
41/// ## Contents
42///
43/// + [Pure component phase equilibria](#pure-component-phase-equilibria)
44/// + [Bubble and dew point calculations](#bubble-and-dew-point-calculations)
45/// + [Flash calculations](#flash-calculations)
46/// + [Heteroazeotropes](#heteroazeotropes)
47/// + [Utility functions](#utility-functions)
48#[derive(Debug, Clone)]
49pub struct PhaseEquilibrium<E, const P: usize, N: Dim = Dyn, D: DualNum<f64> + Copy = f64>
50where
51    DefaultAllocator: Allocator<N>,
52{
53    pub states: [State<E, N, D>; P],
54    pub phase_fractions: [D; P],
55    total_moles: Option<Moles<D>>,
56}
57
58impl<E: Residual<N>, N: Dim, const P: usize> fmt::Display for PhaseEquilibrium<E, P, N>
59where
60    DefaultAllocator: Allocator<N>,
61{
62    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
63        for (i, s) in self.states.iter().enumerate() {
64            writeln!(f, "phase {i}: {s}")?;
65        }
66        Ok(())
67    }
68}
69
70impl<E: Residual, const P: usize> PhaseEquilibrium<E, P> {
71    pub fn _repr_markdown_(&self) -> String {
72        if self.states[0].eos.components() == 1 {
73            let mut res = "||temperature|density|\n|-|-|-|\n".to_string();
74            for (i, s) in self.states.iter().enumerate() {
75                writeln!(
76                    res,
77                    "|phase {}|{:.5}|{:.5}|",
78                    i + 1,
79                    s.temperature,
80                    s.density
81                )
82                .unwrap();
83            }
84            res
85        } else {
86            let mut res = "||temperature|density|molefracs|\n|-|-|-|-|\n".to_string();
87            for (i, s) in self.states.iter().enumerate() {
88                writeln!(
89                    res,
90                    "|phase {}|{:.5}|{:.5}|{:.5?}|",
91                    i + 1,
92                    s.temperature,
93                    s.density,
94                    s.molefracs.as_slice()
95                )
96                .unwrap();
97            }
98            res
99        }
100    }
101}
102
103impl<E: Residual<N, D>, N: Dim, D: DualNum<f64> + Copy> PhaseEquilibrium<E, 2, N, D>
104where
105    DefaultAllocator: Allocator<N>,
106{
107    pub fn vapor(&self) -> &State<E, N, D> {
108        &self.states[0]
109    }
110
111    pub fn liquid(&self) -> &State<E, N, D> {
112        &self.states[1]
113    }
114
115    pub fn vapor_phase_fraction(&self) -> D {
116        self.phase_fractions[0]
117    }
118}
119
120impl<E> PhaseEquilibrium<E, 3> {
121    pub fn vapor(&self) -> &State<E> {
122        &self.states[0]
123    }
124
125    pub fn liquid1(&self) -> &State<E> {
126        &self.states[1]
127    }
128
129    pub fn liquid2(&self) -> &State<E> {
130        &self.states[2]
131    }
132}
133
134impl<E: Residual<N, D>, N: Dim, D: DualNum<f64> + Copy> PhaseEquilibrium<E, 2, N, D>
135where
136    DefaultAllocator: Allocator<N>,
137{
138    pub fn single_phase(state: State<E, N, D>) -> Self {
139        let total_moles = state.total_moles;
140        Self::with_vapor_phase_fraction(state.clone(), state, D::from(1.0), total_moles)
141    }
142
143    pub fn two_phase(vapor: State<E, N, D>, liquid: State<E, N, D>) -> Self {
144        let (beta, total_moles) =
145            if let (Some(nv), Some(nl)) = (vapor.total_moles, liquid.total_moles) {
146                (nv.convert_into(nl + nv), Some(nl + nv))
147            } else {
148                (D::from(1.0), None)
149            };
150        Self::with_vapor_phase_fraction(vapor, liquid, beta, total_moles)
151    }
152
153    pub fn with_vapor_phase_fraction(
154        vapor: State<E, N, D>,
155        liquid: State<E, N, D>,
156        vapor_phase_fraction: D,
157        total_moles: Option<Moles<D>>,
158    ) -> Self {
159        Self {
160            states: [vapor, liquid],
161            phase_fractions: [vapor_phase_fraction, -vapor_phase_fraction + 1.0],
162            total_moles,
163        }
164    }
165}
166
167impl<E: Residual<N, D>, N: Dim, D: DualNum<f64> + Copy> PhaseEquilibrium<E, 3, N, D>
168where
169    DefaultAllocator: Allocator<N>,
170{
171    pub fn new(vapor: State<E, N, D>, liquid1: State<E, N, D>, liquid2: State<E, N, D>) -> Self {
172        Self {
173            states: [vapor, liquid1, liquid2],
174            phase_fractions: [D::from(1.0), D::from(0.0), D::from(0.0)],
175            total_moles: None,
176        }
177    }
178}
179
180impl<E: Residual<N, D>, N: Gradients, const P: usize, D: DualNum<f64> + Copy>
181    PhaseEquilibrium<E, P, N, D>
182where
183    DefaultAllocator: Allocator<N>,
184{
185    pub fn total_moles(&self) -> FeosResult<Moles<D>> {
186        self.total_moles.ok_or(FeosError::IntensiveState)
187    }
188}
189
190impl<E: Total<N, D>, N: Gradients, const P: usize, D: DualNum<f64> + Copy>
191    PhaseEquilibrium<E, P, N, D>
192where
193    DefaultAllocator: Allocator<N>,
194{
195    pub fn molar_enthalpy(&self) -> MolarEnergy<D> {
196        self.states
197            .iter()
198            .zip(&self.phase_fractions)
199            .map(|(s, x)| s.molar_enthalpy(Tot) * Dimensionless::new(x))
200            .reduce(|a, b| a + b)
201            .unwrap()
202    }
203
204    pub fn enthalpy(&self) -> FeosResult<Energy<D>> {
205        Ok(self.total_moles()? * self.molar_enthalpy())
206    }
207
208    pub fn molar_entropy(&self) -> MolarEntropy<D> {
209        self.states
210            .iter()
211            .zip(&self.phase_fractions)
212            .map(|(s, x)| s.molar_entropy(Tot) * Dimensionless::new(x))
213            .reduce(|a, b| a + b)
214            .unwrap()
215    }
216
217    pub fn entropy(&self) -> FeosResult<Entropy<D>> {
218        Ok(self.total_moles()? * self.molar_entropy())
219    }
220}
221
222const TRIVIAL_REL_DEVIATION: f64 = 1e-5;
223
224/// # Utility functions
225impl<E: Residual<N>, N: Dim> PhaseEquilibrium<E, 2, N>
226where
227    DefaultAllocator: Allocator<N>,
228{
229    /// Check if the two states form a trivial solution
230    pub fn is_trivial_solution(state1: &State<E, N>, state2: &State<E, N>) -> bool {
231        let rho1 = state1.molefracs.clone() * state1.density.into_reduced();
232        let rho2 = state2.molefracs.clone() * state2.density.into_reduced();
233
234        rho1.into_iter()
235            .zip(&rho2)
236            .fold(0.0, |acc, (rho1, rho2)| (rho2 / rho1 - 1.0).abs().max(acc))
237            < TRIVIAL_REL_DEVIATION
238    }
239}
feos_core/phase_equilibria/mod.rs

feos_core/phase_equilibria/
mod.rs
