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

feos_core/phase_equilibria/
phase_diagram_binary.rs
