feos_core/phase_equilibria/
phase_diagram_binary.rs

1use super::bubble_dew::TemperatureOrPressure;
2use super::{PhaseDiagram, PhaseEquilibrium};
3use crate::errors::{FeosError, FeosResult};
4use crate::state::{Contributions, DensityInitialization::Vapor, State, StateBuilder};
5use crate::{ReferenceSystem, Residual, SolverOptions, Subset};
6use nalgebra::{DVector, dvector, matrix, stack, vector};
7use ndarray::{Array1, s};
8use num_dual::linalg::LU;
9use quantity::{Density, Moles, Pressure, RGAS, Temperature};
10
11const DEFAULT_POINTS: usize = 51;
12
13impl<E: Residual + Subset> PhaseDiagram<E, 2> {
14    /// Create a new binary phase diagram exhibiting a
15    /// vapor/liquid equilibrium.
16    ///
17    /// If a heteroazeotrope occurs and the composition of the liquid
18    /// phases are known, they can be passed as `x_lle` to avoid
19    /// the calculation of unstable branches.
20    pub fn binary_vle<TP: TemperatureOrPressure>(
21        eos: &E,
22        temperature_or_pressure: TP,
23        npoints: Option<usize>,
24        x_lle: Option<(f64, f64)>,
25        bubble_dew_options: (SolverOptions, SolverOptions),
26    ) -> FeosResult<Self> {
27        let npoints = npoints.unwrap_or(DEFAULT_POINTS);
28
29        // calculate boiling temperature/vapor pressure of pure components
30        let vle_sat = PhaseEquilibrium::vle_pure_comps(eos, temperature_or_pressure);
31        let vle_sat = [vle_sat[1].clone(), vle_sat[0].clone()];
32
33        // Only calculate up to specified compositions
34        if let Some(x_lle) = x_lle {
35            let [states1, states2] = Self::calculate_vlle(
36                eos,
37                temperature_or_pressure,
38                npoints,
39                x_lle,
40                vle_sat,
41                bubble_dew_options,
42            )?;
43
44            let states = states1
45                .into_iter()
46                .chain(states2.into_iter().rev())
47                .collect();
48            return Ok(Self { states });
49        }
50
51        // use dew point when calculating a supercritical tx diagram
52        let bubble = temperature_or_pressure.temperature().is_some();
53
54        // look for supercritical components
55        let (x_lim, vle_lim, bubble) = match vle_sat {
56            [None, None] => return Err(FeosError::SuperCritical),
57            [Some(vle2), None] => {
58                let cp = State::critical_point_binary(
59                    eos,
60                    temperature_or_pressure,
61                    None,
62                    None,
63                    None,
64                    SolverOptions::default(),
65                )?;
66                let cp_vle = PhaseEquilibrium::from_states(cp.clone(), cp.clone());
67                ([0.0, cp.molefracs[0]], (vle2, cp_vle), bubble)
68            }
69            [None, Some(vle1)] => {
70                let cp = State::critical_point_binary(
71                    eos,
72                    temperature_or_pressure,
73                    None,
74                    None,
75                    None,
76                    SolverOptions::default(),
77                )?;
78                let cp_vle = PhaseEquilibrium::from_states(cp.clone(), cp.clone());
79                ([1.0, cp.molefracs[0]], (vle1, cp_vle), bubble)
80            }
81            [Some(vle2), Some(vle1)] => ([0.0, 1.0], (vle2, vle1), true),
82        };
83
84        let mut states = iterate_vle(
85            eos,
86            temperature_or_pressure,
87            &x_lim,
88            vle_lim.0,
89            Some(vle_lim.1),
90            npoints,
91            bubble,
92            bubble_dew_options,
93        );
94        if !bubble {
95            states = states.into_iter().rev().collect();
96        }
97        let states = check_for_vlle(temperature_or_pressure, states, npoints, bubble_dew_options);
98        Ok(Self { states })
99    }
100
101    fn calculate_vlle<TP: TemperatureOrPressure>(
102        eos: &E,
103        tp: TP,
104        npoints: usize,
105        x_lle: (f64, f64),
106        vle_sat: [Option<PhaseEquilibrium<E, 2>>; 2],
107        bubble_dew_options: (SolverOptions, SolverOptions),
108    ) -> FeosResult<[Vec<PhaseEquilibrium<E, 2>>; 2]> {
109        match vle_sat {
110            [Some(vle2), Some(vle1)] => {
111                let states1 = iterate_vle(
112                    eos,
113                    tp,
114                    &[0.0, x_lle.0],
115                    vle2,
116                    None,
117                    npoints / 2,
118                    true,
119                    bubble_dew_options,
120                );
121                let states2 = iterate_vle(
122                    eos,
123                    tp,
124                    &[1.0, x_lle.1],
125                    vle1,
126                    None,
127                    npoints - npoints / 2,
128                    true,
129                    bubble_dew_options,
130                );
131                Ok([states1, states2])
132            }
133            _ => Err(FeosError::SuperCritical),
134        }
135    }
136
137    /// Create a new phase diagram using Tp flash calculations.
138    ///
139    /// The usual use case for this function is the calculation of
140    /// liquid-liquid phase diagrams, but it can be used for vapor-
141    /// liquid diagrams as well, as long as the feed composition is
142    /// in a two phase region.
143    pub fn lle<TP: TemperatureOrPressure>(
144        eos: &E,
145        temperature_or_pressure: TP,
146        feed: &Moles<DVector<f64>>,
147        min_tp: TP::Other,
148        max_tp: TP::Other,
149        npoints: Option<usize>,
150    ) -> FeosResult<Self> {
151        let npoints = npoints.unwrap_or(DEFAULT_POINTS);
152        let mut states = Vec::with_capacity(npoints);
153
154        let (t_vec, p_vec) = temperature_or_pressure.linspace(min_tp, max_tp, npoints);
155        let mut vle = None;
156        for i in 0..npoints {
157            let (t, p) = (t_vec.get(i), p_vec.get(i));
158            vle = PhaseEquilibrium::tp_flash(
159                eos,
160                t,
161                p,
162                feed,
163                vle.as_ref(),
164                SolverOptions::default(),
165                None,
166            )
167            .ok();
168            if let Some(vle) = vle.as_ref() {
169                states.push(vle.clone());
170            }
171        }
172        Ok(Self { states })
173    }
174}
175
176#[expect(clippy::too_many_arguments)]
177fn iterate_vle<E: Residual + Subset, TP: TemperatureOrPressure>(
178    eos: &E,
179    tp: TP,
180    x_lim: &[f64],
181    vle_0: PhaseEquilibrium<E, 2>,
182    vle_1: Option<PhaseEquilibrium<E, 2>>,
183    npoints: usize,
184    bubble: bool,
185    bubble_dew_options: (SolverOptions, SolverOptions),
186) -> Vec<PhaseEquilibrium<E, 2>> {
187    let mut vle_vec = Vec::with_capacity(npoints);
188
189    let x = Array1::linspace(x_lim[0], x_lim[1], npoints);
190    let x = if vle_1.is_some() {
191        x.slice(s![1..-1])
192    } else {
193        x.slice(s![1..])
194    };
195
196    let tp_0 = Some(TP::from_state(vle_0.vapor()));
197    let mut tp_old = tp_0;
198    let mut y_old = None;
199    vle_vec.push(vle_0);
200    for xi in x {
201        let vle = PhaseEquilibrium::bubble_dew_point(
202            eos,
203            tp,
204            &dvector![*xi, 1.0 - xi],
205            tp_old,
206            y_old.as_ref(),
207            bubble,
208            bubble_dew_options,
209        );
210
211        if let Ok(vle) = vle {
212            y_old = Some(if bubble {
213                vle.vapor().molefracs.clone()
214            } else {
215                vle.liquid().molefracs.clone()
216            });
217            tp_old = Some(TP::from_state(vle.vapor()));
218            vle_vec.push(vle.clone());
219        } else {
220            y_old = None;
221            tp_old = tp_0;
222        }
223    }
224    if let Some(vle_1) = vle_1 {
225        vle_vec.push(vle_1);
226    }
227
228    vle_vec
229}
230
231fn check_for_vlle<E: Residual + Subset, TP: TemperatureOrPressure>(
232    tp: TP,
233    states: Vec<PhaseEquilibrium<E, 2>>,
234    npoints: usize,
235    bubble_dew_options: (SolverOptions, SolverOptions),
236) -> Vec<PhaseEquilibrium<E, 2>> {
237    let n = states.len();
238    let p: Vec<_> = states
239        .iter()
240        .map(|s| s.vapor().pressure(Contributions::Total))
241        .collect();
242    let t: Vec<_> = states.iter().map(|s| s.vapor().temperature).collect();
243    let x: Vec<_> = states.iter().map(|s| s.liquid().molefracs[0]).collect();
244    let y: Vec<_> = states.iter().map(|s| s.vapor().molefracs[0]).collect();
245
246    // Determine if the dew line intersects with itself
247    if let Some(t) = tp.temperature()
248        && p[1] > p[0]
249        && p[n - 2] > p[n - 1]
250    {
251        let [mut i, mut j] = [0, n - 1];
252        while i != j {
253            if p[i] > p[j] {
254                j -= 1;
255            } else {
256                i += 1
257            }
258            if y[j] < y[i] {
259                // intersection found!
260                let (xj, yj, pj) = if j == n - 2 {
261                    // Use Henry constant of component 2
262                    let k_inf = (states[n - 1].liquid().ln_phi() - states[n - 1].vapor().ln_phi())
263                        .map(f64::exp)[1];
264                    (
265                        [1.0, 1.0 - 1.0 / k_inf],
266                        [1.0, 0.0],
267                        [p[n - 1], p[n - 1] * (2.0 - 1.0 / k_inf)],
268                    )
269                } else {
270                    // or interpolate linearly
271                    ([x[j + 1], x[j]], [y[j + 1], y[j]], [p[j + 1], p[j]])
272                };
273                let (xi, yi, pi) = if i == 1 {
274                    // Use Henry constant of component 1
275                    let k_inf =
276                        (states[0].liquid().ln_phi() - states[0].vapor().ln_phi()).map(f64::exp)[0];
277                    (
278                        [0.0, 1.0 / k_inf],
279                        [0.0, 1.0],
280                        [p[0], p[0] * (2.0 - 1.0 / k_inf)],
281                    )
282                } else {
283                    // or interpolate linearly
284                    ([x[i - 1], x[i]], [y[i - 1], y[i]], [p[i - 1], p[i]])
285                };
286                // calculate intersection
287                let a = matrix![yi[1] - yi[0], yj[0] - yj[1];
288                                (pi[1] - pi[0]).into_reduced(), (pj[0] - pj[1]).into_reduced()];
289                let b = vector![yj[0] - yi[0], (pj[0] - pi[0]).into_reduced()];
290                let [[r, s]] = LU::new(a).unwrap().solve(&b).data.0;
291                let (xi, xj, p) = (
292                    xi[0] + r * (xi[1] - xi[0]),
293                    xj[0] + s * (xj[1] - xj[0]),
294                    pi[0] + r * (pi[1] - pi[0]),
295                );
296                let Ok(vlle) = PhaseEquilibrium::heteroazeotrope(
297                    &states[0].liquid().eos,
298                    t,
299                    (xi, xj),
300                    Some(p),
301                    Default::default(),
302                    bubble_dew_options,
303                ) else {
304                    return states;
305                };
306                let x_hetero = (vlle.liquid1().molefracs[0], vlle.liquid2().molefracs[0]);
307                return PhaseDiagram::binary_vle(
308                    &states[0].liquid().eos,
309                    tp,
310                    Some(npoints),
311                    Some(x_hetero),
312                    bubble_dew_options,
313                )
314                .map_or(states, |dia| dia.states);
315            }
316        }
317    } else if let Some(p) = tp.pressure()
318        && t[1] < t[0]
319        && t[n - 2] < t[n - 1]
320    {
321        let [mut i, mut j] = [0, n - 1];
322        while i != j {
323            if t[i] < t[j] {
324                j -= 1;
325            } else {
326                i += 1
327            }
328            if y[j] < y[i] {
329                // intersection found!
330                let (xj, yj, tj) = if j == n - 2 {
331                    // Use Henry constant of component 2
332                    let vle = &states[n - 1];
333                    let k_inf = (vle.liquid().ln_phi() - vle.vapor().ln_phi()).map(f64::exp)[1];
334                    let dh = vle.vapor().residual_molar_enthalpy()
335                        - vle.liquid().residual_molar_enthalpy();
336                    let dv = 1.0 / vle.vapor().density - 1.0 / vle.liquid().density;
337                    let pdv_dh = (p * dv).convert_into(dh);
338                    (
339                        [1.0, 1.0 - 1.0 / k_inf],
340                        [1.0, 0.0],
341                        [t[n - 1], t[n - 1] * (1.0 - (k_inf - 1.0) / k_inf * pdv_dh)],
342                    )
343                } else {
344                    // or interpolate linearly
345                    ([x[j + 1], x[j]], [y[j + 1], y[j]], [t[j + 1], t[j]])
346                };
347                let (xi, yi, ti) = if i == 1 {
348                    // Use Henry constant of component 1
349                    let vle = &states[0];
350                    let k_inf = (vle.liquid().ln_phi() - vle.vapor().ln_phi()).map(f64::exp)[0];
351                    let dh = vle.vapor().residual_molar_enthalpy()
352                        - vle.liquid().residual_molar_enthalpy();
353                    let dv = 1.0 / vle.vapor().density - 1.0 / vle.liquid().density;
354                    let pdv_dh = (p * dv).convert_into(dh);
355                    (
356                        [0.0, 1.0 / k_inf],
357                        [0.0, 1.0],
358                        [t[0], t[0] * (1.0 - (k_inf - 1.0) / k_inf * pdv_dh)],
359                    )
360                } else {
361                    // or interpolate linearly
362                    ([x[i - 1], x[i]], [y[i - 1], y[i]], [t[i - 1], t[i]])
363                };
364                // calculate intersection
365                let a = matrix![yi[1] - yi[0], yj[0] - yj[1];
366                                (ti[1] - ti[0]).into_reduced(), (tj[0] - tj[1]).into_reduced()];
367                let b = vector![yj[0] - yi[0], (tj[0] - ti[0]).into_reduced()];
368                let [[r, s]] = LU::new(a).unwrap().solve(&b).data.0;
369                let (xi, xj, t) = (
370                    xi[0] + r * (xi[1] - xi[0]),
371                    xj[0] + s * (xj[1] - xj[0]),
372                    ti[0] + r * (ti[1] - ti[0]),
373                );
374                let Ok(vlle) = PhaseEquilibrium::heteroazeotrope(
375                    &states[0].liquid().eos,
376                    p,
377                    (xi, xj),
378                    Some(t),
379                    Default::default(),
380                    bubble_dew_options,
381                ) else {
382                    return states;
383                };
384                let x_hetero = (vlle.liquid1().molefracs[0], vlle.liquid2().molefracs[0]);
385                return PhaseDiagram::binary_vle(
386                    &states[0].liquid().eos,
387                    tp,
388                    Some(npoints),
389                    Some(x_hetero),
390                    bubble_dew_options,
391                )
392                .map_or(states, |dia| dia.states);
393            }
394        }
395    }
396    states
397}
398
399/// Phase diagram (Txy or pxy) for a system with heteroazeotropic phase behavior.
400pub struct PhaseDiagramHetero<E> {
401    pub vle1: PhaseDiagram<E, 2>,
402    pub vle2: PhaseDiagram<E, 2>,
403    pub lle: Option<PhaseDiagram<E, 2>>,
404}
405
406impl<E: Residual + Subset> PhaseDiagram<E, 2> {
407    /// Create a new binary phase diagram exhibiting a
408    /// vapor/liquid/liquid equilibrium.
409    ///
410    /// The `x_lle` parameter is used as initial values for the calculation
411    /// of the heteroazeotrope.
412    #[expect(clippy::too_many_arguments)]
413    pub fn binary_vlle<TP: TemperatureOrPressure>(
414        eos: &E,
415        temperature_or_pressure: TP,
416        x_lle: (f64, f64),
417        tp_lim_lle: Option<TP::Other>,
418        tp_init_vlle: Option<TP::Other>,
419        npoints_vle: Option<usize>,
420        npoints_lle: Option<usize>,
421        bubble_dew_options: (SolverOptions, SolverOptions),
422    ) -> FeosResult<PhaseDiagramHetero<E>> {
423        let npoints_vle = npoints_vle.unwrap_or(DEFAULT_POINTS);
424
425        // calculate pure components
426        let vle_sat = PhaseEquilibrium::vle_pure_comps(eos, temperature_or_pressure);
427        let vle_sat = [vle_sat[1].clone(), vle_sat[0].clone()];
428
429        // calculate heteroazeotrope
430        let vlle = PhaseEquilibrium::heteroazeotrope(
431            eos,
432            temperature_or_pressure,
433            x_lle,
434            tp_init_vlle,
435            SolverOptions::default(),
436            bubble_dew_options,
437        )?;
438        let x_hetero = (vlle.liquid1().molefracs[0], vlle.liquid2().molefracs[0]);
439
440        // calculate vapor liquid equilibria
441        let [dia1, dia2] = PhaseDiagram::calculate_vlle(
442            eos,
443            temperature_or_pressure,
444            npoints_vle,
445            x_hetero,
446            vle_sat,
447            bubble_dew_options,
448        )?;
449
450        // calculate liquid liquid equilibrium
451        let lle = tp_lim_lle
452            .map(|tp_lim| {
453                let tp_hetero = TP::from_state(vlle.vapor());
454                let x_feed = 0.5 * (x_hetero.0 + x_hetero.1);
455                let feed = Moles::from_reduced(dvector![x_feed, 1.0 - x_feed]);
456                PhaseDiagram::lle(
457                    eos,
458                    temperature_or_pressure,
459                    &feed,
460                    tp_lim,
461                    tp_hetero,
462                    npoints_lle,
463                )
464            })
465            .transpose()?;
466
467        Ok(PhaseDiagramHetero {
468            vle1: PhaseDiagram::new(dia1),
469            vle2: PhaseDiagram::new(dia2),
470            lle,
471        })
472    }
473}
474
475impl<E: Clone> PhaseDiagramHetero<E> {
476    pub fn vle(&self) -> PhaseDiagram<E, 2> {
477        PhaseDiagram::new(
478            self.vle1
479                .states
480                .iter()
481                .chain(self.vle2.states.iter().rev())
482                .cloned()
483                .collect(),
484        )
485    }
486}
487
488const MAX_ITER_HETERO: usize = 50;
489const TOL_HETERO: f64 = 1e-8;
490
491/// # Heteroazeotropes
492impl<E: Residual> PhaseEquilibrium<E, 3> {
493    /// Calculate a heteroazeotrope (three phase equilbrium) for a binary
494    /// system and given temperature or pressure.
495    pub fn heteroazeotrope<TP: TemperatureOrPressure>(
496        eos: &E,
497        temperature_or_pressure: TP,
498        x_init: (f64, f64),
499        tp_init: Option<TP::Other>,
500        options: SolverOptions,
501        bubble_dew_options: (SolverOptions, SolverOptions),
502    ) -> FeosResult<Self> {
503        let (temperature, pressure, iterate_p) =
504            temperature_or_pressure.temperature_pressure(tp_init);
505        if iterate_p {
506            PhaseEquilibrium::heteroazeotrope_t(
507                eos,
508                temperature.unwrap(),
509                x_init,
510                pressure,
511                options,
512                bubble_dew_options,
513            )
514        } else {
515            PhaseEquilibrium::heteroazeotrope_p(
516                eos,
517                pressure.unwrap(),
518                x_init,
519                temperature,
520                options,
521                bubble_dew_options,
522            )
523        }
524    }
525
526    /// Calculate a heteroazeotrope (three phase equilbrium) for a binary
527    /// system and given temperature.
528    #[expect(clippy::toplevel_ref_arg)]
529    fn heteroazeotrope_t(
530        eos: &E,
531        temperature: Temperature,
532        x_init: (f64, f64),
533        p_init: Option<Pressure>,
534        options: SolverOptions,
535        bubble_dew_options: (SolverOptions, SolverOptions),
536    ) -> FeosResult<Self> {
537        // calculate initial values using bubble point
538        let x1 = dvector![x_init.0, 1.0 - x_init.0];
539        let x2 = dvector![x_init.1, 1.0 - x_init.1];
540        let vle1 = PhaseEquilibrium::bubble_point(
541            eos,
542            temperature,
543            &x1,
544            p_init,
545            None,
546            bubble_dew_options,
547        )?;
548        let vle2 = PhaseEquilibrium::bubble_point(
549            eos,
550            temperature,
551            &x2,
552            p_init,
553            None,
554            bubble_dew_options,
555        )?;
556        let mut l1 = vle1.liquid().clone();
557        let mut l2 = vle2.liquid().clone();
558        let p0 = (vle1.vapor().pressure(Contributions::Total)
559            + vle2.vapor().pressure(Contributions::Total))
560            * 0.5;
561        let nv0 = (&vle1.vapor().moles + &vle2.vapor().moles) * 0.5;
562        let mut v = State::new_npt(eos, temperature, p0, &nv0, Some(Vapor))?;
563
564        for _ in 0..options.max_iter.unwrap_or(MAX_ITER_HETERO) {
565            // calculate properties
566            let dmu_drho_l1 = (l1.dmu_dni(Contributions::Total) * l1.volume).to_reduced();
567            let dmu_drho_l2 = (l2.dmu_dni(Contributions::Total) * l2.volume).to_reduced();
568            let dmu_drho_v = (v.dmu_dni(Contributions::Total) * v.volume).to_reduced();
569            let dp_drho_l1 = (l1.dp_dni(Contributions::Total) * l1.volume)
570                .to_reduced()
571                .transpose();
572            let dp_drho_l2 = (l2.dp_dni(Contributions::Total) * l2.volume)
573                .to_reduced()
574                .transpose();
575            let dp_drho_v = (v.dp_dni(Contributions::Total) * v.volume)
576                .to_reduced()
577                .transpose();
578            let mu_l1_res = l1.residual_chemical_potential().to_reduced();
579            let mu_l2_res = l2.residual_chemical_potential().to_reduced();
580            let mu_v_res = v.residual_chemical_potential().to_reduced();
581            let p_l1 = l1.pressure(Contributions::Total).to_reduced();
582            let p_l2 = l2.pressure(Contributions::Total).to_reduced();
583            let p_v = v.pressure(Contributions::Total).to_reduced();
584
585            // calculate residual
586            let delta_l1v_mu_ig = (RGAS * v.temperature).to_reduced()
587                * (l1
588                    .partial_density
589                    .to_reduced()
590                    .component_div(&v.partial_density.to_reduced()))
591                .map(f64::ln);
592            let delta_l2v_mu_ig = (RGAS * v.temperature).to_reduced()
593                * (l2
594                    .partial_density
595                    .to_reduced()
596                    .component_div(&v.partial_density.to_reduced()))
597                .map(f64::ln);
598            let res = stack![
599                mu_l1_res - &mu_v_res + delta_l1v_mu_ig;
600                mu_l2_res - &mu_v_res + delta_l2v_mu_ig;
601                vector![p_l1 - p_v];
602                vector![p_l2 - p_v]
603            ];
604
605            // check for convergence
606            if res.norm() < options.tol.unwrap_or(TOL_HETERO) {
607                return Ok(Self([v, l1, l2]));
608            }
609
610            // calculate Jacobian
611            let jacobian = stack![
612                dmu_drho_l1, 0          , -&dmu_drho_v;
613                0          , dmu_drho_l2, -dmu_drho_v;
614                dp_drho_l1 , 0          , -&dp_drho_v;
615                0          , dp_drho_l2 , -dp_drho_v
616            ];
617
618            // calculate Newton step
619            let dx = LU::new(jacobian)?.solve(&res);
620
621            // apply Newton step
622            let rho_l1 =
623                &l1.partial_density - &Density::from_reduced(dx.rows_range(0..2).into_owned());
624            let rho_l2 =
625                &l2.partial_density - &Density::from_reduced(dx.rows_range(2..4).into_owned());
626            let rho_v =
627                &v.partial_density - &Density::from_reduced(dx.rows_range(4..6).into_owned());
628
629            // check for negative densities
630            for i in 0..2 {
631                if rho_l1.get(i).is_sign_negative()
632                    || rho_l2.get(i).is_sign_negative()
633                    || rho_v.get(i).is_sign_negative()
634                {
635                    return Err(FeosError::IterationFailed(String::from(
636                        "PhaseEquilibrium::heteroazeotrope_t",
637                    )));
638                }
639            }
640
641            // update states
642            l1 = StateBuilder::new(eos)
643                .temperature(temperature)
644                .partial_density(&rho_l1)
645                .build()?;
646            l2 = StateBuilder::new(eos)
647                .temperature(temperature)
648                .partial_density(&rho_l2)
649                .build()?;
650            v = StateBuilder::new(eos)
651                .temperature(temperature)
652                .partial_density(&rho_v)
653                .build()?;
654        }
655        Err(FeosError::NotConverged(String::from(
656            "PhaseEquilibrium::heteroazeotrope_t",
657        )))
658    }
659
660    /// Calculate a heteroazeotrope (three phase equilbrium) for a binary
661    /// system and given pressure.
662    #[expect(clippy::toplevel_ref_arg)]
663    fn heteroazeotrope_p(
664        eos: &E,
665        pressure: Pressure,
666        x_init: (f64, f64),
667        t_init: Option<Temperature>,
668        options: SolverOptions,
669        bubble_dew_options: (SolverOptions, SolverOptions),
670    ) -> FeosResult<Self> {
671        let p = pressure.to_reduced();
672
673        // calculate initial values using bubble point
674        let x1 = dvector![x_init.0, 1.0 - x_init.0];
675        let x2 = dvector![x_init.1, 1.0 - x_init.1];
676        let vle1 =
677            PhaseEquilibrium::bubble_point(eos, pressure, &x1, t_init, None, bubble_dew_options)?;
678        let vle2 =
679            PhaseEquilibrium::bubble_point(eos, pressure, &x2, t_init, None, bubble_dew_options)?;
680        let mut l1 = vle1.liquid().clone();
681        let mut l2 = vle2.liquid().clone();
682        let t0 = (vle1.vapor().temperature + vle2.vapor().temperature) * 0.5;
683        let nv0 = (&vle1.vapor().moles + &vle2.vapor().moles) * 0.5;
684        let mut v = State::new_npt(eos, t0, pressure, &nv0, Some(Vapor))?;
685
686        for _ in 0..options.max_iter.unwrap_or(MAX_ITER_HETERO) {
687            // calculate properties
688            let dmu_drho_l1 = (l1.dmu_dni(Contributions::Total) * l1.volume).to_reduced();
689            let dmu_drho_l2 = (l2.dmu_dni(Contributions::Total) * l2.volume).to_reduced();
690            let dmu_drho_v = (v.dmu_dni(Contributions::Total) * v.volume).to_reduced();
691            let dmu_res_dt_l1 = (l1.dmu_res_dt()).to_reduced();
692            let dmu_res_dt_l2 = (l2.dmu_res_dt()).to_reduced();
693            let dmu_res_dt_v = (v.dmu_res_dt()).to_reduced();
694            let dp_drho_l1 = (l1.dp_dni(Contributions::Total) * l1.volume)
695                .to_reduced()
696                .transpose();
697            let dp_drho_l2 = (l2.dp_dni(Contributions::Total) * l2.volume)
698                .to_reduced()
699                .transpose();
700            let dp_drho_v = (v.dp_dni(Contributions::Total) * v.volume)
701                .to_reduced()
702                .transpose();
703            let dp_dt_l1 = (l1.dp_dt(Contributions::Total)).to_reduced();
704            let dp_dt_l2 = (l2.dp_dt(Contributions::Total)).to_reduced();
705            let dp_dt_v = (v.dp_dt(Contributions::Total)).to_reduced();
706            let mu_l1_res = l1.residual_chemical_potential().to_reduced();
707            let mu_l2_res = l2.residual_chemical_potential().to_reduced();
708            let mu_v_res = v.residual_chemical_potential().to_reduced();
709            let p_l1 = l1.pressure(Contributions::Total).to_reduced();
710            let p_l2 = l2.pressure(Contributions::Total).to_reduced();
711            let p_v = v.pressure(Contributions::Total).to_reduced();
712
713            // calculate residual
714            let delta_l1v_dmu_ig_dt = l1
715                .partial_density
716                .to_reduced()
717                .component_div(&v.partial_density.to_reduced())
718                .map(f64::ln);
719            let delta_l2v_dmu_ig_dt = l2
720                .partial_density
721                .to_reduced()
722                .component_div(&v.partial_density.to_reduced())
723                .map(f64::ln);
724            let delta_l1v_mu_ig = (RGAS * v.temperature).to_reduced() * &delta_l1v_dmu_ig_dt;
725            let delta_l2v_mu_ig = (RGAS * v.temperature).to_reduced() * &delta_l2v_dmu_ig_dt;
726            let res = stack![
727                mu_l1_res - &mu_v_res + delta_l1v_mu_ig;
728                mu_l2_res - &mu_v_res + delta_l2v_mu_ig;
729                vector![p_l1 - p];
730                vector![p_l2 - p];
731                vector![p_v - p]
732            ];
733
734            // check for convergence
735            if res.norm() < options.tol.unwrap_or(TOL_HETERO) {
736                return Ok(Self([v, l1, l2]));
737            }
738
739            let jacobian = stack![
740                dmu_drho_l1, 0, -&dmu_drho_v, dmu_res_dt_l1 - &dmu_res_dt_v + delta_l1v_dmu_ig_dt;
741                0, dmu_drho_l2, -dmu_drho_v, dmu_res_dt_l2 - &dmu_res_dt_v + delta_l2v_dmu_ig_dt;
742                dp_drho_l1, 0, 0, vector![dp_dt_l1];
743                0, dp_drho_l2, 0, vector![dp_dt_l2];
744                0, 0, dp_drho_v, vector![dp_dt_v]
745            ];
746
747            // calculate Newton step
748            let dx = LU::new(jacobian)?.solve(&res);
749
750            // apply Newton step
751            let rho_l1 =
752                l1.partial_density - Density::from_reduced(dx.rows_range(0..2).into_owned());
753            let rho_l2 =
754                l2.partial_density - Density::from_reduced(dx.rows_range(2..4).into_owned());
755            let rho_v = v.partial_density - Density::from_reduced(dx.rows_range(4..6).into_owned());
756            let t = v.temperature - Temperature::from_reduced(dx[6]);
757
758            // check for negative densities and temperatures
759            for i in 0..2 {
760                if rho_l1.get(i).is_sign_negative()
761                    || rho_l2.get(i).is_sign_negative()
762                    || rho_v.get(i).is_sign_negative()
763                    || t.is_sign_negative()
764                {
765                    return Err(FeosError::IterationFailed(String::from(
766                        "PhaseEquilibrium::heteroazeotrope_p",
767                    )));
768                }
769            }
770
771            // update states
772            l1 = StateBuilder::new(eos)
773                .temperature(t)
774                .partial_density(&rho_l1)
775                .build()?;
776            l2 = StateBuilder::new(eos)
777                .temperature(t)
778                .partial_density(&rho_l2)
779                .build()?;
780            v = StateBuilder::new(eos)
781                .temperature(t)
782                .partial_density(&rho_v)
783                .build()?;
784        }
785        Err(FeosError::NotConverged(String::from(
786            "PhaseEquilibrium::heteroazeotrope_p",
787        )))
788    }
789}
feos_core/phase_equilibria/phase_diagram_binary.rs

feos_core/phase_equilibria/
phase_diagram_binary.rs
