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, 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        Ok(Self { states })
98    }
99
100    #[expect(clippy::type_complexity)]
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>>, Vec<PhaseEquilibrium<E, 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
231/// Phase diagram (Txy or pxy) for a system with heteroazeotropic phase behavior.
232pub struct PhaseDiagramHetero<E> {
233    pub vle1: PhaseDiagram<E, 2>,
234    pub vle2: PhaseDiagram<E, 2>,
235    pub lle: Option<PhaseDiagram<E, 2>>,
236}
237
238impl<E: Residual + Subset> PhaseDiagram<E, 2> {
239    /// Create a new binary phase diagram exhibiting a
240    /// vapor/liquid/liquid equilibrium.
241    ///
242    /// The `x_lle` parameter is used as initial values for the calculation
243    /// of the heteroazeotrope.
244    #[expect(clippy::too_many_arguments)]
245    pub fn binary_vlle<TP: TemperatureOrPressure>(
246        eos: &E,
247        temperature_or_pressure: TP,
248        x_lle: (f64, f64),
249        tp_lim_lle: Option<TP::Other>,
250        tp_init_vlle: Option<TP::Other>,
251        npoints_vle: Option<usize>,
252        npoints_lle: Option<usize>,
253        bubble_dew_options: (SolverOptions, SolverOptions),
254    ) -> FeosResult<PhaseDiagramHetero<E>> {
255        let npoints_vle = npoints_vle.unwrap_or(DEFAULT_POINTS);
256
257        // calculate pure components
258        let vle_sat = PhaseEquilibrium::vle_pure_comps(eos, temperature_or_pressure);
259        let vle_sat = [vle_sat[1].clone(), vle_sat[0].clone()];
260
261        // calculate heteroazeotrope
262        let vlle = PhaseEquilibrium::heteroazeotrope(
263            eos,
264            temperature_or_pressure,
265            x_lle,
266            tp_init_vlle,
267            SolverOptions::default(),
268            bubble_dew_options,
269        )?;
270        let x_hetero = (vlle.liquid1().molefracs[0], vlle.liquid2().molefracs[0]);
271
272        // calculate vapor liquid equilibria
273        let (dia1, dia2) = PhaseDiagram::calculate_vlle(
274            eos,
275            temperature_or_pressure,
276            npoints_vle,
277            x_hetero,
278            vle_sat,
279            bubble_dew_options,
280        )?;
281
282        // calculate liquid liquid equilibrium
283        let lle = tp_lim_lle
284            .map(|tp_lim| {
285                let tp_hetero = TP::from_state(vlle.vapor());
286                let x_feed = 0.5 * (x_hetero.0 + x_hetero.1);
287                let feed = Moles::from_reduced(dvector![x_feed, 1.0 - x_feed]);
288                PhaseDiagram::lle(
289                    eos,
290                    temperature_or_pressure,
291                    &feed,
292                    tp_lim,
293                    tp_hetero,
294                    npoints_lle,
295                )
296            })
297            .transpose()?;
298
299        Ok(PhaseDiagramHetero {
300            vle1: PhaseDiagram::new(dia1),
301            vle2: PhaseDiagram::new(dia2),
302            lle,
303        })
304    }
305}
306
307impl<E: Clone> PhaseDiagramHetero<E> {
308    pub fn vle(&self) -> PhaseDiagram<E, 2> {
309        PhaseDiagram::new(
310            self.vle1
311                .states
312                .iter()
313                .chain(self.vle2.states.iter().rev())
314                .cloned()
315                .collect(),
316        )
317    }
318}
319
320const MAX_ITER_HETERO: usize = 50;
321const TOL_HETERO: f64 = 1e-8;
322
323/// # Heteroazeotropes
324impl<E: Residual> PhaseEquilibrium<E, 3> {
325    /// Calculate a heteroazeotrope (three phase equilbrium) for a binary
326    /// system and given temperature or pressure.
327    pub fn heteroazeotrope<TP: TemperatureOrPressure>(
328        eos: &E,
329        temperature_or_pressure: TP,
330        x_init: (f64, f64),
331        tp_init: Option<TP::Other>,
332        options: SolverOptions,
333        bubble_dew_options: (SolverOptions, SolverOptions),
334    ) -> FeosResult<Self> {
335        let (temperature, pressure, iterate_p) =
336            temperature_or_pressure.temperature_pressure(tp_init);
337        // let tp_init = tp_init.map(|tp_init| temperature_or_pressure.temperature_pressure(tp_init));
338        if iterate_p {
339            PhaseEquilibrium::heteroazeotrope_t(
340                eos,
341                temperature.unwrap(),
342                x_init,
343                pressure,
344                options,
345                bubble_dew_options,
346            )
347        } else {
348            PhaseEquilibrium::heteroazeotrope_p(
349                eos,
350                pressure.unwrap(),
351                x_init,
352                temperature,
353                options,
354                bubble_dew_options,
355            )
356        }
357    }
358
359    /// Calculate a heteroazeotrope (three phase equilbrium) for a binary
360    /// system and given temperature.
361    #[expect(clippy::toplevel_ref_arg)]
362    fn heteroazeotrope_t(
363        eos: &E,
364        temperature: Temperature,
365        x_init: (f64, f64),
366        p_init: Option<Pressure>,
367        options: SolverOptions,
368        bubble_dew_options: (SolverOptions, SolverOptions),
369    ) -> FeosResult<Self> {
370        // calculate initial values using bubble point
371        let x1 = dvector![x_init.0, 1.0 - x_init.0];
372        let x2 = dvector![x_init.1, 1.0 - x_init.1];
373        let vle1 = PhaseEquilibrium::bubble_point(
374            eos,
375            temperature,
376            &x1,
377            p_init,
378            None,
379            bubble_dew_options,
380        )?;
381        let vle2 = PhaseEquilibrium::bubble_point(
382            eos,
383            temperature,
384            &x2,
385            p_init,
386            None,
387            bubble_dew_options,
388        )?;
389        let mut l1 = vle1.liquid().clone();
390        let mut l2 = vle2.liquid().clone();
391        let p0 = (vle1.vapor().pressure(Contributions::Total)
392            + vle2.vapor().pressure(Contributions::Total))
393            * 0.5;
394        let nv0 = (&vle1.vapor().moles + &vle2.vapor().moles) * 0.5;
395        let mut v = State::new_npt(eos, temperature, p0, &nv0, Some(Vapor))?;
396
397        for _ in 0..options.max_iter.unwrap_or(MAX_ITER_HETERO) {
398            // calculate properties
399            let dmu_drho_l1 = (l1.dmu_dni(Contributions::Total) * l1.volume).to_reduced();
400            let dmu_drho_l2 = (l2.dmu_dni(Contributions::Total) * l2.volume).to_reduced();
401            let dmu_drho_v = (v.dmu_dni(Contributions::Total) * v.volume).to_reduced();
402            let dp_drho_l1 = (l1.dp_dni(Contributions::Total) * l1.volume)
403                .to_reduced()
404                .transpose();
405            let dp_drho_l2 = (l2.dp_dni(Contributions::Total) * l2.volume)
406                .to_reduced()
407                .transpose();
408            let dp_drho_v = (v.dp_dni(Contributions::Total) * v.volume)
409                .to_reduced()
410                .transpose();
411            let mu_l1_res = l1.residual_chemical_potential().to_reduced();
412            let mu_l2_res = l2.residual_chemical_potential().to_reduced();
413            let mu_v_res = v.residual_chemical_potential().to_reduced();
414            let p_l1 = l1.pressure(Contributions::Total).to_reduced();
415            let p_l2 = l2.pressure(Contributions::Total).to_reduced();
416            let p_v = v.pressure(Contributions::Total).to_reduced();
417
418            // calculate residual
419            let delta_l1v_mu_ig = (RGAS * v.temperature).to_reduced()
420                * (l1
421                    .partial_density
422                    .to_reduced()
423                    .component_div(&v.partial_density.to_reduced()))
424                .map(f64::ln);
425            let delta_l2v_mu_ig = (RGAS * v.temperature).to_reduced()
426                * (l2
427                    .partial_density
428                    .to_reduced()
429                    .component_div(&v.partial_density.to_reduced()))
430                .map(f64::ln);
431            let res = stack![
432                mu_l1_res - &mu_v_res + delta_l1v_mu_ig;
433                mu_l2_res - &mu_v_res + delta_l2v_mu_ig;
434                vector![p_l1 - p_v];
435                vector![p_l2 - p_v]
436            ];
437
438            // check for convergence
439            if res.norm() < options.tol.unwrap_or(TOL_HETERO) {
440                return Ok(Self([v, l1, l2]));
441            }
442
443            // calculate Jacobian
444            let jacobian = stack![
445                    dmu_drho_l1, 0, 0;
446                    0, dmu_drho_l2, -dmu_drho_v;
447                    dp_drho_l1, 0, -&dp_drho_v;
448                    0, dp_drho_l2, -dp_drho_v
449            ];
450
451            // calculate Newton step
452            let dx = LU::new(jacobian)?.solve(&res);
453
454            // apply Newton step
455            let rho_l1 =
456                &l1.partial_density - &Density::from_reduced(dx.rows_range(0..2).into_owned());
457            let rho_l2 =
458                &l2.partial_density - &Density::from_reduced(dx.rows_range(2..4).into_owned());
459            let rho_v =
460                &v.partial_density - &Density::from_reduced(dx.rows_range(4..6).into_owned());
461
462            // check for negative densities
463            for i in 0..2 {
464                if rho_l1.get(i).is_sign_negative()
465                    || rho_l2.get(i).is_sign_negative()
466                    || rho_v.get(i).is_sign_negative()
467                {
468                    return Err(FeosError::IterationFailed(String::from(
469                        "PhaseEquilibrium::heteroazeotrope_t",
470                    )));
471                }
472            }
473
474            // update states
475            l1 = StateBuilder::new(eos)
476                .temperature(temperature)
477                .partial_density(&rho_l1)
478                .build()?;
479            l2 = StateBuilder::new(eos)
480                .temperature(temperature)
481                .partial_density(&rho_l2)
482                .build()?;
483            v = StateBuilder::new(eos)
484                .temperature(temperature)
485                .partial_density(&rho_v)
486                .build()?;
487        }
488        Err(FeosError::NotConverged(String::from(
489            "PhaseEquilibrium::heteroazeotrope_t",
490        )))
491    }
492
493    /// Calculate a heteroazeotrope (three phase equilbrium) for a binary
494    /// system and given pressure.
495    #[expect(clippy::toplevel_ref_arg)]
496    fn heteroazeotrope_p(
497        eos: &E,
498        pressure: Pressure,
499        x_init: (f64, f64),
500        t_init: Option<Temperature>,
501        options: SolverOptions,
502        bubble_dew_options: (SolverOptions, SolverOptions),
503    ) -> FeosResult<Self> {
504        let p = pressure.to_reduced();
505
506        // calculate initial values using bubble point
507        let x1 = dvector![x_init.0, 1.0 - x_init.0];
508        let x2 = dvector![x_init.1, 1.0 - x_init.1];
509        let vle1 =
510            PhaseEquilibrium::bubble_point(eos, pressure, &x1, t_init, None, bubble_dew_options)?;
511        let vle2 =
512            PhaseEquilibrium::bubble_point(eos, pressure, &x2, t_init, None, bubble_dew_options)?;
513        let mut l1 = vle1.liquid().clone();
514        let mut l2 = vle2.liquid().clone();
515        let t0 = (vle1.vapor().temperature + vle2.vapor().temperature) * 0.5;
516        let nv0 = (&vle1.vapor().moles + &vle2.vapor().moles) * 0.5;
517        let mut v = State::new_npt(eos, t0, pressure, &nv0, Some(Vapor))?;
518
519        for _ in 0..options.max_iter.unwrap_or(MAX_ITER_HETERO) {
520            // calculate properties
521            let dmu_drho_l1 = (l1.dmu_dni(Contributions::Total) * l1.volume).to_reduced();
522            let dmu_drho_l2 = (l2.dmu_dni(Contributions::Total) * l2.volume).to_reduced();
523            let dmu_drho_v = (v.dmu_dni(Contributions::Total) * v.volume).to_reduced();
524            let dmu_res_dt_l1 = (l1.dmu_res_dt()).to_reduced();
525            let dmu_res_dt_l2 = (l2.dmu_res_dt()).to_reduced();
526            let dmu_res_dt_v = (v.dmu_res_dt()).to_reduced();
527            let dp_drho_l1 = (l1.dp_dni(Contributions::Total) * l1.volume)
528                .to_reduced()
529                .transpose();
530            let dp_drho_l2 = (l2.dp_dni(Contributions::Total) * l2.volume)
531                .to_reduced()
532                .transpose();
533            let dp_drho_v = (v.dp_dni(Contributions::Total) * v.volume)
534                .to_reduced()
535                .transpose();
536            let dp_dt_l1 = (l1.dp_dt(Contributions::Total)).to_reduced();
537            let dp_dt_l2 = (l2.dp_dt(Contributions::Total)).to_reduced();
538            let dp_dt_v = (v.dp_dt(Contributions::Total)).to_reduced();
539            let mu_l1_res = l1.residual_chemical_potential().to_reduced();
540            let mu_l2_res = l2.residual_chemical_potential().to_reduced();
541            let mu_v_res = v.residual_chemical_potential().to_reduced();
542            let p_l1 = l1.pressure(Contributions::Total).to_reduced();
543            let p_l2 = l2.pressure(Contributions::Total).to_reduced();
544            let p_v = v.pressure(Contributions::Total).to_reduced();
545
546            // calculate residual
547            let delta_l1v_dmu_ig_dt = l1
548                .partial_density
549                .to_reduced()
550                .component_div(&v.partial_density.to_reduced())
551                .map(f64::ln);
552            let delta_l2v_dmu_ig_dt = l2
553                .partial_density
554                .to_reduced()
555                .component_div(&v.partial_density.to_reduced())
556                .map(f64::ln);
557            let delta_l1v_mu_ig = (RGAS * v.temperature).to_reduced() * &delta_l1v_dmu_ig_dt;
558            let delta_l2v_mu_ig = (RGAS * v.temperature).to_reduced() * &delta_l2v_dmu_ig_dt;
559            let res = stack![
560                mu_l1_res - &mu_v_res + delta_l1v_mu_ig;
561                mu_l2_res - &mu_v_res + delta_l2v_mu_ig;
562                vector![p_l1 - p];
563                vector![p_l2 - p];
564                vector![p_v - p]
565            ];
566
567            // check for convergence
568            if res.norm() < options.tol.unwrap_or(TOL_HETERO) {
569                return Ok(Self([v, l1, l2]));
570            }
571
572            let jacobian = stack![
573                dmu_drho_l1, 0, -&dmu_drho_v, dmu_res_dt_l1 - &dmu_res_dt_v + delta_l1v_dmu_ig_dt;
574                0, dmu_drho_l2, -dmu_drho_v, dmu_res_dt_l2 - &dmu_res_dt_v + delta_l2v_dmu_ig_dt;
575                dp_drho_l1, 0, 0, vector![dp_dt_l1];
576                0, dp_drho_l2, 0, vector![dp_dt_l2];
577                0, 0, dp_drho_v, vector![dp_dt_v]
578            ];
579
580            // calculate Newton step
581            let dx = LU::new(jacobian)?.solve(&res);
582
583            // apply Newton step
584            let rho_l1 =
585                l1.partial_density - Density::from_reduced(dx.rows_range(0..2).into_owned());
586            let rho_l2 =
587                l2.partial_density - Density::from_reduced(dx.rows_range(2..4).into_owned());
588            let rho_v = v.partial_density - Density::from_reduced(dx.rows_range(4..6).into_owned());
589            let t = v.temperature - Temperature::from_reduced(dx[6]);
590
591            // check for negative densities and temperatures
592            for i in 0..2 {
593                if rho_l1.get(i).is_sign_negative()
594                    || rho_l2.get(i).is_sign_negative()
595                    || rho_v.get(i).is_sign_negative()
596                    || t.is_sign_negative()
597                {
598                    return Err(FeosError::IterationFailed(String::from(
599                        "PhaseEquilibrium::heteroazeotrope_p",
600                    )));
601                }
602            }
603
604            // update states
605            l1 = StateBuilder::new(eos)
606                .temperature(t)
607                .partial_density(&rho_l1)
608                .build()?;
609            l2 = StateBuilder::new(eos)
610                .temperature(t)
611                .partial_density(&rho_l2)
612                .build()?;
613            v = StateBuilder::new(eos)
614                .temperature(t)
615                .partial_density(&rho_v)
616                .build()?;
617        }
618        Err(FeosError::NotConverged(String::from(
619            "PhaseEquilibrium::heteroazeotrope_p",
620        )))
621    }
622}
feos_core/phase_equilibria/phase_diagram_binary.rs

feos_core/phase_equilibria/
phase_diagram_binary.rs
