feos_ad/core/
phase_equilibria.rs

1use super::{HelmholtzEnergyWrapper, ParametersAD, ResidualHelmholtzEnergy, StateAD};
2use feos_core::{Contributions, EosResult, PhaseEquilibrium, ReferenceSystem};
3use nalgebra::SVector;
4use ndarray::{arr1, Array};
5use num_dual::{linalg::LU, DualNum};
6use quantity::{Dimensionless, Moles, Pressure, Temperature};
7
8impl<'a, R: ResidualHelmholtzEnergy<N>, D: DualNum<f64> + Copy, const N: usize>
9    StateAD<'a, R, D, N>
10{
11    /// Perform a Tp-flash calculation. Returns the [PhaseEquilibriumAD] and the vapor fraction.
12    pub fn tp_flash(&self) -> EosResult<(PhaseEquilibriumAD<'a, R, D, N>, Dimensionless<D>)> {
13        let pressure = self.pressure();
14        let feed = Moles::from_reduced(arr1(&self.molefracs.data.0[0].map(|x| x.re())));
15        let vle = PhaseEquilibrium::tp_flash(
16            &self.eos.eos,
17            self.temperature.re(),
18            pressure.re(),
19            &feed,
20            None,
21            Default::default(),
22            None,
23        )?;
24        let rho_l = vle.liquid().partial_density.to_reduced();
25        let mut rho_l = SVector::from_fn(|i, _| D::from(rho_l[i]));
26        let rho_v = vle.vapor().partial_density.to_reduced();
27        let mut rho_v = SVector::from_fn(|i, _| D::from(rho_v[i]));
28        let mut v_l = D::from(vle.liquid().volume.to_reduced());
29        let mut v_v = D::from(vle.vapor().volume.to_reduced());
30        let t = self.reduced_temperature;
31        let p = pressure.into_reduced();
32
33        for _ in 0..D::NDERIV {
34            let (p_l, mu_res_l, dp_l, dmu_l) = R::dmu_drho(&self.eos.parameters, t, &rho_l);
35            let (p_v, mu_res_v, dp_v, dmu_v) = R::dmu_drho(&self.eos.parameters, t, &rho_v);
36
37            let f = Array::from_shape_fn((2 * N + 2,), |i| {
38                if i < N {
39                    mu_res_l[i] - mu_res_v[i] + (rho_l[i] / rho_v[i]).ln() * t
40                } else if i < 2 * N {
41                    rho_l[i - N] * v_l + rho_v[i - N] * v_v - self.molefracs[i - N]
42                } else if i == 2 * N {
43                    p_l - p
44                } else if i == 2 * N + 1 {
45                    p_v - p
46                } else {
47                    unreachable!()
48                }
49            });
50            let jac = Array::from_shape_fn((2 * N + 2, 2 * N + 2), |(i, j)| {
51                if i < N {
52                    if j < N {
53                        dmu_l[(i, j)]
54                    } else if j < 2 * N {
55                        -dmu_v[(i, j - N)]
56                    } else {
57                        D::zero()
58                    }
59                } else if i < 2 * N {
60                    if j + N == i {
61                        v_l
62                    } else if j == i {
63                        v_v
64                    } else if j == 2 * N {
65                        rho_l[i - N]
66                    } else if j == 2 * N + 1 {
67                        rho_v[i - N]
68                    } else {
69                        D::zero()
70                    }
71                } else if i == 2 * N && j < N {
72                    dp_l[j]
73                } else if i == 2 * N + 1 && N <= j && j < 2 * N {
74                    dp_v[j - N]
75                } else {
76                    D::zero()
77                }
78            });
79
80            let dx = LU::new(jac).unwrap().solve(&(-f));
81            let drho_l = SVector::from_fn(|i, _| dx[i]);
82            let drho_v = SVector::from_fn(|i, _| dx[i + N]);
83            let dv_l = dx[2 * N];
84            let dv_v = dx[2 * N + 1];
85
86            rho_l += drho_l;
87            rho_v += drho_v;
88            v_l += dv_l;
89            v_v += dv_v;
90        }
91        let molar_volume_l = rho_l.sum().recip();
92        let molar_volume_v = rho_v.sum().recip();
93        let molefracs_l = rho_l * molar_volume_l;
94        let molefracs_v = rho_v * molar_volume_v;
95        Ok((
96            PhaseEquilibriumAD {
97                liquid: StateAD::new(self.eos, t, molar_volume_l, molefracs_l),
98                vapor: StateAD::new(self.eos, t, molar_volume_v, molefracs_v),
99            },
100            Dimensionless::from_reduced(v_v / molar_volume_v / self.molefracs.sum()),
101        ))
102    }
103}
104
105/// An equilibrium state consisting of a vapor and a liquid phase.
106pub struct PhaseEquilibriumAD<'a, E: ParametersAD, D: DualNum<f64> + Copy, const N: usize> {
107    pub liquid: StateAD<'a, E, D, N>,
108    pub vapor: StateAD<'a, E, D, N>,
109}
110
111impl<'a, R: ResidualHelmholtzEnergy<1>, D: DualNum<f64> + Copy> PhaseEquilibriumAD<'a, R, D, 1> {
112    /// Calculate a phase equilibrium of a pure component for a given temperature.
113    /// Returns the phase equilibrium and the vapor pressure.
114    pub fn new_t(
115        eos: &'a HelmholtzEnergyWrapper<R, D, 1>,
116        temperature: Temperature<D>,
117    ) -> EosResult<(Self, Pressure<D>)> {
118        let vle = PhaseEquilibrium::pure(&eos.eos, temperature.re(), None, Default::default())?;
119        let mut density1 = D::from(vle.liquid().density.to_reduced());
120        let mut density2 = D::from(vle.vapor().density.to_reduced());
121        let molefracs = SVector::from([D::one()]);
122        let t = temperature.into_reduced();
123        let mut p = D::from(vle.vapor().pressure(Contributions::Total).to_reduced());
124        for _ in 0..D::NDERIV {
125            let (f1, p1, dp_drho1) = R::dp_drho(&eos.parameters, t, density1.recip(), &molefracs);
126            let (f2, p2, dp_drho2) = R::dp_drho(&eos.parameters, t, density2.recip(), &molefracs);
127            p = -(density2 * f1 - density1 * f2
128                + density1 * density2 * t * (density1 / density2).ln())
129                / (density2 - density1);
130            density1 -= (p1 - p) / dp_drho1;
131            density2 -= (p2 - p) / dp_drho2;
132        }
133        Ok((
134            Self {
135                liquid: StateAD::new(eos, t, density1.recip(), molefracs),
136                vapor: StateAD::new(eos, t, density2.recip(), molefracs),
137            },
138            Pressure::from_reduced(p),
139        ))
140    }
141}
142
143impl<'a, R: ResidualHelmholtzEnergy<N>, D: DualNum<f64> + Copy, const N: usize>
144    PhaseEquilibriumAD<'a, R, D, N>
145{
146    /// Calculate a bubble point of a mixture for a given temperature.
147    /// Returns the phase equilibrium and the bubble point pressure.
148    pub fn bubble_point(
149        eos: &'a HelmholtzEnergyWrapper<R, D, N>,
150        temperature: Temperature<D>,
151        liquid_molefracs: SVector<D, N>,
152    ) -> EosResult<(Self, Pressure<D>)> {
153        let x = arr1(liquid_molefracs.map(|x| x.re()).as_slice());
154        let vle = PhaseEquilibrium::bubble_point(
155            &eos.eos,
156            temperature.re(),
157            &x,
158            None,
159            None,
160            Default::default(),
161        )?;
162        let rho_v = vle.vapor().partial_density.to_reduced();
163        let (liquid, vapor, pressure) = Self::bubble_dew_point(
164            eos,
165            temperature,
166            liquid_molefracs,
167            vle.liquid().pressure(Contributions::Total).to_reduced(),
168            vle.liquid().density.to_reduced(),
169            SVector::from_fn(|i, _| rho_v[i]),
170        )?;
171        Ok((Self { liquid, vapor }, pressure))
172    }
173
174    /// Calculate a dew point of a mixture for a given temperature.
175    /// Returns the phase equilibrium and the dew point pressure.
176    pub fn dew_point(
177        eos: &'a HelmholtzEnergyWrapper<R, D, N>,
178        temperature: Temperature<D>,
179        vapor_molefracs: SVector<D, N>,
180    ) -> EosResult<(Self, Pressure<D>)> {
181        let y = arr1(vapor_molefracs.map(|y| y.re()).as_slice());
182        let vle = PhaseEquilibrium::dew_point(
183            &eos.eos,
184            temperature.re(),
185            &y,
186            None,
187            None,
188            Default::default(),
189        )?;
190        let rho_l = vle.liquid().partial_density.to_reduced();
191        let (vapor, liquid, pressure) = Self::bubble_dew_point(
192            eos,
193            temperature,
194            vapor_molefracs,
195            vle.vapor().pressure(Contributions::Total).to_reduced(),
196            vle.vapor().density.to_reduced(),
197            SVector::from_fn(|i, _| rho_l[i]),
198        )?;
199        Ok((Self { liquid, vapor }, pressure))
200    }
201
202    #[expect(clippy::type_complexity)]
203    fn bubble_dew_point(
204        eos: &'a HelmholtzEnergyWrapper<R, D, N>,
205        temperature: Temperature<D>,
206        molefracs: SVector<D, N>,
207        pressure: f64,
208        density: f64,
209        partial_density_other_phase: SVector<f64, N>,
210    ) -> EosResult<(StateAD<'a, R, D, N>, StateAD<'a, R, D, N>, Pressure<D>)> {
211        let mut rho = SVector::from_fn(|i, _| D::from(partial_density_other_phase[i]));
212        let mut v = D::from(density.recip());
213        let t = temperature.into_reduced();
214        let mut p = D::from(pressure);
215        for _ in 0..D::NDERIV {
216            let (p_1, mu_res_1, dp_1, dmu_1) = R::dmu_drho(&eos.parameters, t, &rho);
217            let (p_2, mu_res_2, dp_2, dmu_2) = R::dmu_dv(&eos.parameters, t, v, &molefracs);
218
219            let f = Array::from_shape_fn((N + 2,), |i| {
220                if i == N {
221                    p_1 - p
222                } else if i == N + 1 {
223                    p_2 - p
224                } else {
225                    mu_res_1[i] - mu_res_2[i] + (rho[i] * v / molefracs[i]).ln() * t
226                }
227            });
228            let jac = Array::from_shape_fn((N + 2, N + 2), |(i, j)| {
229                if i < N && j < N {
230                    dmu_1[(i, j)]
231                } else if i < N && j == N {
232                    -dmu_2[i]
233                } else if i == N && j < N {
234                    dp_1[j]
235                } else if i == N + 1 && j == N {
236                    dp_2
237                } else if i >= N && j == N + 1 {
238                    -D::one()
239                } else {
240                    D::zero()
241                }
242            });
243
244            let dx = LU::new(jac).unwrap().solve(&(-f));
245            let drho = SVector::from_fn(|i, _| dx[i]);
246            let dv = dx[N];
247            let dp = dx[N + 1];
248
249            rho += drho;
250            v += dv;
251            p += dp;
252        }
253        let v_o = rho.sum().recip();
254        let molefracs_other_phase = rho * v_o;
255        Ok((
256            StateAD::new(eos, t, v, molefracs),
257            StateAD::new(eos, t, v_o, molefracs_other_phase),
258            Pressure::from_reduced(p),
259        ))
260    }
261}
262
263#[cfg(test)]
264mod test {
265    use super::*;
266    use crate::eos::pcsaft::test::pcsaft;
267    use crate::eos::{GcPcSaft, GcPcSaftParameters, Joback};
268    use crate::EquationOfStateAD;
269    use approx::assert_relative_eq;
270    use feos_core::{Contributions, DensityInitialization, EosResult, PhaseEquilibrium};
271    use num_dual::{Dual, Dual64};
272    use quantity::KELVIN;
273    use std::collections::HashMap;
274
275    #[test]
276    fn test_phase_equilibrium() -> EosResult<()> {
277        let (pcsaft, eos) = pcsaft()?;
278        let pcsaft = pcsaft.wrap();
279        let temperature = 250.0 * KELVIN;
280        let (vle, p) = PhaseEquilibriumAD::new_t(&pcsaft, temperature)?;
281        let rho_l = 1.0 / vle.liquid.molar_volume;
282        let rho_v = 1.0 / vle.vapor.molar_volume;
283        let vle_feos = PhaseEquilibrium::pure(&eos, temperature, None, Default::default())?;
284        let p_feos = vle_feos.vapor().pressure(Contributions::Total);
285        println!("{:.5} {:.5} {:.5}", rho_l, rho_v, p);
286        println!(
287            "{:.5} {:.5} {:.5}",
288            vle_feos.liquid().density,
289            vle_feos.vapor().density,
290            p_feos
291        );
292        println!("{} {}", vle.liquid.pressure(), vle.vapor.pressure());
293        assert_relative_eq!(rho_l, vle_feos.liquid().density, max_relative = 1e-10);
294        assert_relative_eq!(rho_v, vle_feos.vapor().density, max_relative = 1e-10);
295        assert_relative_eq!(p, p_feos, max_relative = 1e-10);
296        Ok(())
297    }
298
299    #[test]
300    fn test_phase_equilibrium_derivative() -> EosResult<()> {
301        let (pcsaft, eos) = pcsaft()?;
302        let eos_ad = pcsaft.wrap().derivatives(pcsaft.params());
303        let t: Temperature<Dual64> = Temperature::from_reduced(Dual::from(250.0).derivative());
304        let (vle, p) = PhaseEquilibriumAD::new_t(&eos_ad, t)?;
305        let rho_l = vle.liquid.reduced_molar_volume.recip();
306        let rho_v = vle.vapor.reduced_molar_volume.recip();
307        let p = p.into_reduced();
308        let vle_feos = PhaseEquilibrium::pure(&eos, 250.0 * KELVIN, None, Default::default())?;
309        let h = 1e-5 * KELVIN;
310        let vle_feos_h =
311            PhaseEquilibrium::pure(&eos, 250.0 * KELVIN + h, None, Default::default())?;
312        let drho_l = ((vle_feos_h.liquid().density - vle_feos.liquid().density) / h).to_reduced();
313        let drho_v = ((vle_feos_h.vapor().density - vle_feos.vapor().density) / h).to_reduced();
314        let dp = ((vle_feos_h.vapor().pressure(Contributions::Total)
315            - vle_feos.vapor().pressure(Contributions::Total))
316            / h)
317            .to_reduced();
318        println!("{:11.5e} {:11.5e} {:11.5e}", rho_l.eps, rho_v.eps, p.eps);
319        println!("{:11.5e} {:11.5e} {:11.5e}", drho_l, drho_v, dp,);
320        println!(
321            "{} {}",
322            vle.liquid.pressure().into_reduced(),
323            vle.vapor.pressure().into_reduced()
324        );
325        assert_relative_eq!(rho_l.eps, drho_l, max_relative = 1e-5);
326        assert_relative_eq!(rho_v.eps, drho_v, max_relative = 1e-5);
327        assert_relative_eq!(p.eps, dp, max_relative = 1e-5);
328        Ok(())
329    }
330
331    fn acetone_pentane_parameters() -> GcPcSaftParameters<f64, 2> {
332        let mut groups1 = HashMap::new();
333        groups1.insert("CH3", 2.0);
334        groups1.insert(">C=O", 1.0);
335        let mut bonds1 = HashMap::new();
336        bonds1.insert(["CH3", ">C=O"], 2.0);
337        let mut groups2 = HashMap::new();
338        groups2.insert("CH3", 2.0);
339        groups2.insert("CH2", 3.0);
340        let mut bonds2 = HashMap::new();
341        bonds2.insert(["CH3", "CH2"], 2.0);
342        bonds2.insert(["CH2", "CH2"], 2.0);
343
344        GcPcSaftParameters::from_groups([&groups1, &groups2], [&bonds1, &bonds2])
345    }
346
347    fn acetone_groups() -> HashMap<&'static str, f64> {
348        let mut groups = HashMap::new();
349        groups.insert("CH3", 2.0);
350        groups.insert(">C=O", 1.0);
351        groups
352    }
353
354    fn pentane_groups() -> HashMap<&'static str, f64> {
355        let mut groups = HashMap::new();
356        groups.insert("CH3", 2.0);
357        groups.insert("CH2", 3.0);
358        groups
359    }
360
361    #[test]
362    fn test_dew_point() -> EosResult<()> {
363        let params = GcPcSaft(acetone_pentane_parameters());
364        let joback = [
365            Joback(Joback::from_group_counts(&acetone_groups())),
366            Joback(Joback::from_group_counts(&pentane_groups())),
367        ];
368
369        let mut params_dual = params.params::<Dual64>();
370        params_dual.groups[0].eps = 1.0;
371        let joback_dual = joback.map(|j| j.params());
372
373        let mut params_h = GcPcSaft(acetone_pentane_parameters());
374        let h = 1e-7;
375        params_h.0.groups[0] += h;
376
377        let eos = EquationOfStateAD::new(joback, params).wrap();
378        let eos_h = EquationOfStateAD::new(joback, params_h).wrap();
379        let eos_dual = eos.derivatives((joback_dual, params_dual));
380
381        let temperature = Temperature::from_reduced(Dual::from(300.0));
382        let vapor_molefracs = SVector::from([Dual::from(0.5); 2]);
383        let (vle, p_dew) = PhaseEquilibriumAD::dew_point(
384            &eos,
385            Temperature::from_reduced(300.0),
386            SVector::from([0.5, 0.5]),
387        )?;
388        let (vle_h, p_dew_h) = PhaseEquilibriumAD::dew_point(
389            &eos_h,
390            Temperature::from_reduced(300.0),
391            SVector::from([0.5, 0.5]),
392        )?;
393        let (vle_dual, p_dew_dual) =
394            PhaseEquilibriumAD::dew_point(&eos_dual, temperature, vapor_molefracs)?;
395
396        println!(
397            "{:.6} + {:.6}ε   {:.6} + {:.6}ε   {:.6} + {:.6}ε",
398            vle.vapor.reduced_molar_volume,
399            (vle_h.vapor.reduced_molar_volume - vle.vapor.reduced_molar_volume) / h,
400            vle.liquid.reduced_molar_volume,
401            (vle_h.liquid.reduced_molar_volume - vle.liquid.reduced_molar_volume) / h,
402            p_dew.into_reduced(),
403            (p_dew_h.into_reduced() - p_dew.into_reduced()) / h
404        );
405        println!(
406            "{:.6} + {:.6}ε   {:.6} + {:.6}ε   {:.6} + {:.6}ε",
407            vle_dual.vapor.reduced_molar_volume.re,
408            vle_dual.vapor.reduced_molar_volume.eps,
409            vle_dual.liquid.reduced_molar_volume.re,
410            vle_dual.liquid.reduced_molar_volume.eps,
411            p_dew_dual.into_reduced().re,
412            p_dew_dual.into_reduced().eps,
413        );
414
415        println!(
416            "{:.6} + {:.6}ε   {:.6} + {:.6}ε",
417            vle.liquid.molefracs[0],
418            (vle_h.liquid.molefracs[0] - vle.liquid.molefracs[0]) / h,
419            vle.liquid.molefracs[1],
420            (vle_h.liquid.molefracs[1] - vle.liquid.molefracs[1]) / h,
421        );
422        println!(
423            "{:.6} + {:.6}ε   {:.6} + {:.6}ε",
424            vle_dual.liquid.molefracs[0].re,
425            vle_dual.liquid.molefracs[0].eps,
426            vle_dual.liquid.molefracs[1].re,
427            vle_dual.liquid.molefracs[1].eps
428        );
429
430        let dx = (vle_h.liquid.molefracs[0] - vle.liquid.molefracs[0]) / h;
431        assert_relative_eq!(vle_dual.liquid.molefracs[0].eps, dx, max_relative = 1e-6);
432        let dp = (p_dew_h.into_reduced() - p_dew.into_reduced()) / h;
433        assert_relative_eq!(p_dew_dual.into_reduced().eps, dp, max_relative = 1e-6);
434
435        Ok(())
436    }
437
438    #[test]
439    fn test_bubble_point() -> EosResult<()> {
440        let params = GcPcSaft(acetone_pentane_parameters());
441        let joback = [
442            Joback(Joback::from_group_counts(&acetone_groups())),
443            Joback(Joback::from_group_counts(&pentane_groups())),
444        ];
445
446        let mut params_dual = params.params::<Dual64>();
447        params_dual.groups[0].eps = 1.0;
448        let joback_dual = joback.map(|j| j.params());
449
450        let mut params_h = GcPcSaft(acetone_pentane_parameters());
451        let h = 1e-7;
452        params_h.0.groups[0] += h;
453
454        let eos = EquationOfStateAD::new(joback, params).wrap();
455        let eos_h = EquationOfStateAD::new(joback, params_h).wrap();
456        let eos_dual = eos.derivatives((joback_dual, params_dual));
457
458        let temperature = Temperature::from_reduced(Dual::from(300.0));
459        let liquid_molefracs = SVector::from([Dual::from(0.5); 2]);
460        let (vle, p_bubble) = PhaseEquilibriumAD::bubble_point(
461            &eos,
462            Temperature::from_reduced(300.0),
463            SVector::from([0.5, 0.5]),
464        )?;
465        let (vle_h, p_bubble_h) = PhaseEquilibriumAD::bubble_point(
466            &eos_h,
467            Temperature::from_reduced(300.0),
468            SVector::from([0.5, 0.5]),
469        )?;
470        let (vle_dual, p_bubble_dual) =
471            PhaseEquilibriumAD::bubble_point(&eos_dual, temperature, liquid_molefracs)?;
472
473        println!(
474            "{:.6} + {:.6}ε   {:.6} + {:.6}ε   {:.6} + {:.6}ε",
475            vle.vapor.reduced_molar_volume,
476            (vle_h.vapor.reduced_molar_volume - vle.vapor.reduced_molar_volume) / h,
477            vle.liquid.reduced_molar_volume,
478            (vle_h.liquid.reduced_molar_volume - vle.liquid.reduced_molar_volume) / h,
479            p_bubble.into_reduced(),
480            (p_bubble_h.into_reduced() - p_bubble.into_reduced()) / h
481        );
482        println!(
483            "{:.6} + {:.6}ε   {:.6} + {:.6}ε   {:.6} + {:.6}ε",
484            vle_dual.vapor.reduced_molar_volume.re,
485            vle_dual.vapor.reduced_molar_volume.eps,
486            vle_dual.liquid.reduced_molar_volume.re,
487            vle_dual.liquid.reduced_molar_volume.eps,
488            p_bubble_dual.into_reduced().re,
489            p_bubble_dual.into_reduced().eps,
490        );
491
492        println!(
493            "{:.6} + {:.6}ε   {:.6} + {:.6}ε",
494            vle.vapor.molefracs[0],
495            (vle_h.vapor.molefracs[0] - vle.vapor.molefracs[0]) / h,
496            vle.vapor.molefracs[1],
497            (vle_h.vapor.molefracs[1] - vle.vapor.molefracs[1]) / h,
498        );
499        println!(
500            "{:.6} + {:.6}ε   {:.6} + {:.6}ε",
501            vle_dual.vapor.molefracs[0].re,
502            vle_dual.vapor.molefracs[0].eps,
503            vle_dual.vapor.molefracs[1].re,
504            vle_dual.vapor.molefracs[1].eps
505        );
506
507        let dx = (vle_h.vapor.molefracs[0] - vle.vapor.molefracs[0]) / h;
508        assert_relative_eq!(vle_dual.vapor.molefracs[0].eps, dx, max_relative = 1e-6);
509        let dp = (p_bubble_h.into_reduced() - p_bubble.into_reduced()) / h;
510        assert_relative_eq!(p_bubble_dual.into_reduced().eps, dp, max_relative = 1e-3);
511
512        Ok(())
513    }
514
515    #[test]
516    fn test_tp_flash() -> EosResult<()> {
517        let params = GcPcSaft(acetone_pentane_parameters());
518        let joback = [
519            Joback(Joback::from_group_counts(&acetone_groups())),
520            Joback(Joback::from_group_counts(&pentane_groups())),
521        ];
522
523        let mut params_dual = params.params::<Dual64>();
524        params_dual.groups[0].eps = 1.0;
525        let joback_dual = joback.map(|j| j.params());
526
527        let mut params_h = GcPcSaft(acetone_pentane_parameters());
528        let h = 1e-5;
529        params_h.0.groups[0] += h;
530
531        let eos = EquationOfStateAD::new(joback, params).wrap();
532        let eos_h = EquationOfStateAD::new(joback, params_h).wrap();
533        let eos_dual = eos.derivatives((joback_dual, params_dual));
534
535        let temperature = Temperature::from_reduced(Dual::from(300.0));
536        let pressure = Pressure::from_reduced(Dual::from(0.005));
537        let molefracs = SVector::from([Dual::from(0.5); 2]);
538        let (vle_dual, phi_dual) = StateAD::new_tp(
539            &eos_dual,
540            temperature,
541            pressure,
542            molefracs,
543            DensityInitialization::None,
544        )?
545        .tp_flash()?;
546        let (vle, phi) = StateAD::new_tp(
547            &eos,
548            Temperature::from_reduced(300.0),
549            Pressure::from_reduced(0.005),
550            SVector::from([0.5, 0.5]),
551            DensityInitialization::None,
552        )?
553        .tp_flash()?;
554        let (vle_h, phi_h) = StateAD::new_tp(
555            &eos_h,
556            Temperature::from_reduced(300.0),
557            Pressure::from_reduced(0.005),
558            SVector::from([0.5, 0.5]),
559            DensityInitialization::None,
560        )?
561        .tp_flash()?;
562
563        println!(
564            "{:.6} + {:.6}ε   {:.6} + {:.6}ε   {:.6} + {:.6}ε",
565            vle.vapor.reduced_molar_volume,
566            (vle_h.vapor.reduced_molar_volume - vle.vapor.reduced_molar_volume) / h,
567            vle.liquid.reduced_molar_volume,
568            (vle_h.liquid.reduced_molar_volume - vle.liquid.reduced_molar_volume) / h,
569            phi.to_reduced(),
570            (phi_h - phi).to_reduced() / h
571        );
572        println!(
573            "{:.6} + {:.6}ε   {:.6} + {:.6}ε   {:.6} + {:.6}ε",
574            vle_dual.vapor.reduced_molar_volume.re,
575            vle_dual.vapor.reduced_molar_volume.eps,
576            vle_dual.liquid.reduced_molar_volume.re,
577            vle_dual.liquid.reduced_molar_volume.eps,
578            phi_dual.into_reduced().re,
579            phi_dual.into_reduced().eps
580        );
581
582        println!(
583            "{:.6} + {:.6}ε   {:.6} + {:.6}ε",
584            vle.vapor.molefracs[0],
585            (vle_h.vapor.molefracs[0] - vle.vapor.molefracs[0]) / h,
586            vle.vapor.molefracs[1],
587            (vle_h.vapor.molefracs[1] - vle.vapor.molefracs[1]) / h,
588        );
589        println!(
590            "{:.6} + {:.6}ε   {:.6} + {:.6}ε",
591            vle_dual.vapor.molefracs[0].re,
592            vle_dual.vapor.molefracs[0].eps,
593            vle_dual.vapor.molefracs[1].re,
594            vle_dual.vapor.molefracs[1].eps
595        );
596
597        println!(
598            "{:.6} + {:.6}ε   {:.6} + {:.6}ε",
599            vle.liquid.molefracs[0],
600            (vle_h.liquid.molefracs[0] - vle.liquid.molefracs[0]) / h,
601            vle.liquid.molefracs[1],
602            (vle_h.liquid.molefracs[1] - vle.liquid.molefracs[1]) / h,
603        );
604        println!(
605            "{:.6} + {:.6}ε   {:.6} + {:.6}ε",
606            vle_dual.liquid.molefracs[0].re,
607            vle_dual.liquid.molefracs[0].eps,
608            vle_dual.liquid.molefracs[1].re,
609            vle_dual.liquid.molefracs[1].eps
610        );
611
612        let dx = (vle_h.vapor.molefracs[0] - vle.vapor.molefracs[0]) / h;
613        assert_relative_eq!(vle_dual.vapor.molefracs[0].eps, dx, max_relative = 1e-4);
614        assert_relative_eq!(
615            phi_dual.into_reduced().eps,
616            (phi_h - phi).into_reduced() / h,
617            max_relative = 1e-4
618        );
619
620        Ok(())
621    }
622}
feos_ad/core/phase_equilibria.rs

feos_ad/core/
phase_equilibria.rs
