feos_core/phase_equilibria/
mod.rs

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

feos_core/phase_equilibria/
mod.rs
