feos_core/phase_equilibria/
phase_envelope.rs

1use super::{PhaseDiagram, PhaseEquilibrium};
2use crate::SolverOptions;
3use crate::equation_of_state::Residual;
4use crate::errors::FeosResult;
5use crate::state::{Contributions, State};
6use nalgebra::DVector;
7use quantity::{Pressure, Temperature};
8
9impl<E: Residual> PhaseDiagram<E, 2> {
10    /// Calculate the bubble point line of a mixture with given composition.
11    pub fn bubble_point_line(
12        eos: &E,
13        molefracs: &DVector<f64>,
14        min_temperature: Temperature,
15        npoints: usize,
16        critical_temperature: Option<Temperature>,
17        options: (SolverOptions, SolverOptions),
18    ) -> FeosResult<Self> {
19        let mut states = Vec::with_capacity(npoints);
20
21        let sc = State::critical_point(
22            eos,
23            Some(molefracs),
24            critical_temperature,
25            None,
26            SolverOptions::default(),
27        )?;
28
29        let max_temperature = min_temperature
30            + (sc.temperature - min_temperature) * ((npoints - 2) as f64 / (npoints - 1) as f64);
31        let temperatures = Temperature::linspace(min_temperature, max_temperature, npoints - 1);
32
33        let mut vle: Option<PhaseEquilibrium<E, 2>> = None;
34        for ti in &temperatures {
35            // calculate new liquid point
36            let p_init = vle
37                .as_ref()
38                .map(|vle| vle.vapor().pressure(Contributions::Total));
39            let vapor_molefracs = vle.as_ref().map(|vle| &vle.vapor().molefracs);
40            vle = PhaseEquilibrium::bubble_point(
41                eos,
42                ti,
43                molefracs,
44                p_init,
45                vapor_molefracs,
46                options,
47            )
48            .ok();
49
50            if let Some(vle) = vle.as_ref() {
51                states.push(vle.clone());
52            }
53        }
54        states.push(PhaseEquilibrium::from_states(sc.clone(), sc));
55
56        Ok(PhaseDiagram::new(states))
57    }
58
59    /// Calculate the dew point line of a mixture with given composition.
60    pub fn dew_point_line(
61        eos: &E,
62        molefracs: &DVector<f64>,
63        min_temperature: Temperature,
64        npoints: usize,
65        critical_temperature: Option<Temperature>,
66        options: (SolverOptions, SolverOptions),
67    ) -> FeosResult<Self> {
68        let mut states = Vec::with_capacity(npoints);
69
70        let sc = State::critical_point(
71            eos,
72            Some(molefracs),
73            critical_temperature,
74            None,
75            SolverOptions::default(),
76        )?;
77
78        let n_t = npoints / 2;
79        let max_temperature = min_temperature
80            + (sc.temperature - min_temperature) * ((n_t - 2) as f64 / (n_t - 1) as f64);
81        let temperatures = Temperature::linspace(min_temperature, max_temperature, n_t - 1);
82
83        let mut vle: Option<PhaseEquilibrium<E, 2>> = None;
84        for ti in &temperatures {
85            let p_init = vle
86                .as_ref()
87                .map(|vle| vle.vapor().pressure(Contributions::Total));
88            let liquid_molefracs = vle.as_ref().map(|vle| &vle.liquid().molefracs);
89            vle =
90                PhaseEquilibrium::dew_point(eos, ti, molefracs, p_init, liquid_molefracs, options)
91                    .ok();
92            if let Some(vle) = vle.as_ref() {
93                states.push(vle.clone());
94            }
95        }
96
97        let n_p = npoints - n_t;
98        if vle.is_none() {
99            return Ok(PhaseDiagram::new(states));
100        }
101
102        let min_pressure = vle.as_ref().unwrap().vapor().pressure(Contributions::Total);
103        let p_c = sc.pressure(Contributions::Total);
104        let max_pressure =
105            min_pressure + (p_c - min_pressure) * ((n_p - 2) as f64 / (n_p - 1) as f64);
106        let pressures = Pressure::linspace(min_pressure, max_pressure, n_p);
107
108        for pi in &pressures {
109            let t_init = vle.as_ref().map(|vle| vle.vapor().temperature);
110            let liquid_molefracs = vle.as_ref().map(|vle| &vle.liquid().molefracs);
111            vle =
112                PhaseEquilibrium::dew_point(eos, pi, molefracs, t_init, liquid_molefracs, options)
113                    .ok();
114            if let Some(vle) = vle.as_ref() {
115                states.push(vle.clone());
116            }
117        }
118
119        states.push(PhaseEquilibrium::from_states(sc.clone(), sc));
120
121        Ok(PhaseDiagram::new(states))
122    }
123
124    /// Calculate the spinodal lines for a mixture with fixed composition.
125    pub fn spinodal(
126        eos: &E,
127        molefracs: &DVector<f64>,
128        min_temperature: Temperature,
129        npoints: usize,
130        critical_temperature: Option<Temperature>,
131        options: SolverOptions,
132    ) -> FeosResult<Self> {
133        let mut states = Vec::with_capacity(npoints);
134
135        let sc = State::critical_point(
136            eos,
137            Some(molefracs),
138            critical_temperature,
139            None,
140            SolverOptions::default(),
141        )?;
142
143        let max_temperature = min_temperature
144            + (sc.temperature - min_temperature) * ((npoints - 2) as f64 / (npoints - 1) as f64);
145        let temperatures = Temperature::linspace(min_temperature, max_temperature, npoints - 1);
146
147        for ti in &temperatures {
148            let spinodal = State::spinodal(eos, ti, Some(molefracs), options).ok();
149            if let Some(spinodal) = spinodal {
150                states.push(PhaseEquilibrium(spinodal));
151            }
152        }
153        states.push(PhaseEquilibrium::from_states(sc.clone(), sc));
154
155        Ok(PhaseDiagram::new(states))
156    }
157}
feos_core/phase_equilibria/phase_envelope.rs

feos_core/phase_equilibria/
phase_envelope.rs
