feos_core/phase_equilibria/
phase_envelope.rs

1use super::{PhaseDiagram, PhaseEquilibrium};
2use crate::equation_of_state::Residual;
3use crate::errors::FeosResult;
4use crate::state::{Contributions, State};
5use crate::{Composition, SolverOptions};
6use nalgebra::Dyn;
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<X: Composition<f64, Dyn> + Clone>(
12        eos: &E,
13        composition: X,
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            composition.clone(),
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                composition.clone(),
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::single_phase(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<X: Composition<f64, Dyn> + Clone>(
61        eos: &E,
62        composition: X,
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            composition.clone(),
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 = PhaseEquilibrium::dew_point(
90                eos,
91                ti,
92                composition.clone(),
93                p_init,
94                liquid_molefracs,
95                options,
96            )
97            .ok();
98            if let Some(vle) = vle.as_ref() {
99                states.push(vle.clone());
100            }
101        }
102
103        let n_p = npoints - n_t;
104        if vle.is_none() {
105            return Ok(PhaseDiagram::new(states));
106        }
107
108        let min_pressure = vle.as_ref().unwrap().vapor().pressure(Contributions::Total);
109        let p_c = sc.pressure(Contributions::Total);
110        let max_pressure =
111            min_pressure + (p_c - min_pressure) * ((n_p - 2) as f64 / (n_p - 1) as f64);
112        let pressures = Pressure::linspace(min_pressure, max_pressure, n_p);
113
114        for pi in &pressures {
115            let t_init = vle.as_ref().map(|vle| vle.vapor().temperature);
116            let liquid_molefracs = vle.as_ref().map(|vle| &vle.liquid().molefracs);
117            vle = PhaseEquilibrium::dew_point(
118                eos,
119                pi,
120                composition.clone(),
121                t_init,
122                liquid_molefracs,
123                options,
124            )
125            .ok();
126            if let Some(vle) = vle.as_ref() {
127                states.push(vle.clone());
128            }
129        }
130
131        states.push(PhaseEquilibrium::single_phase(sc));
132
133        Ok(PhaseDiagram::new(states))
134    }
135
136    /// Calculate the spinodal lines for a mixture with fixed composition.
137    pub fn spinodal<X: Composition<f64, Dyn> + Clone>(
138        eos: &E,
139        composition: X,
140        min_temperature: Temperature,
141        npoints: usize,
142        critical_temperature: Option<Temperature>,
143        options: SolverOptions,
144    ) -> FeosResult<Self> {
145        let mut states = Vec::with_capacity(npoints);
146
147        let sc = State::critical_point(
148            eos,
149            composition.clone(),
150            critical_temperature,
151            None,
152            SolverOptions::default(),
153        )?;
154
155        let max_temperature = min_temperature
156            + (sc.temperature - min_temperature) * ((npoints - 2) as f64 / (npoints - 1) as f64);
157        let temperatures = Temperature::linspace(min_temperature, max_temperature, npoints - 1);
158
159        for ti in &temperatures {
160            let spinodal = State::spinodal(eos, ti, composition.clone(), options).ok();
161            if let Some([sp_v, sp_l]) = spinodal {
162                states.push(PhaseEquilibrium::two_phase(sp_v, sp_l));
163            }
164        }
165        states.push(PhaseEquilibrium::single_phase(sc));
166
167        Ok(PhaseDiagram::new(states))
168    }
169}
feos_core/phase_equilibria/phase_envelope.rs

feos_core/phase_equilibria/
phase_envelope.rs
