diff --git a/CHANGELOG.md b/CHANGELOG.md index c0502f7ee..93e16be91 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,9 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). ## [Unreleased] +## Fixed +- Fixed the calculation of critical points and tp-flashes when one or more of the mixture components have zero composition. [#331](https://github.com/feos-org/feos/pull/331) + ### Packaging - Updated `quantity` dependency to 0.13 and removed the `typenum` dependency. [#323](https://github.com/feos-org/feos/pull/323) diff --git a/crates/feos-core/src/phase_equilibria/tp_flash.rs b/crates/feos-core/src/phase_equilibria/tp_flash.rs index 1c9c3072a..baefea80b 100644 --- a/crates/feos-core/src/phase_equilibria/tp_flash.rs +++ b/crates/feos-core/src/phase_equilibria/tp_flash.rs @@ -265,17 +265,13 @@ impl PhaseEquilibrium { // check for convergence *iter += 1; - let mut res_vec = ln_phi_l - ln_phi_v - + self - .liquid() - .molefracs - .component_div(&self.vapor().molefracs) - .map(|i| if i > 0.0 { i.ln() } else { 0.0 }); + let mut res_vec = k.component_mul(&self.liquid().molefracs) - &self.vapor().molefracs; // Set residuum to 0 for non-volatile components if let Some(nvc) = non_volatile_components.as_ref() { nvc.iter().for_each(|&c| res_vec[c] = 0.0); } + let res = res_vec.norm(); log_iter!( verbosity, diff --git a/crates/feos-core/src/state/critical_point.rs b/crates/feos-core/src/state/critical_point.rs index 1f18093e0..e77b95299 100644 --- a/crates/feos-core/src/state/critical_point.rs +++ b/crates/feos-core/src/state/critical_point.rs @@ -10,6 +10,7 @@ use num_dual::{ implicit_derivative_binary, implicit_derivative_vec, jacobian, partial, partial2, third_derivative, }; +use num_traits::Zero; use quantity::{Density, Pressure, Temperature}; const MAX_ITER_CRIT_POINT: usize = 50; @@ -537,7 +538,13 @@ where let n = n + sqrt_z.component_mul(&u).map(Dual3::from_re) * s; let t = Dual3::from_re(temperature); let v = Dual3::from_re(molar_volume); - let ig = n.dot(&n.map(|n| (n / v).ln() - 1.0)); + let ig = n.dot(&n.map(|n| { + if n.is_zero() { + Dual3::zero() + } else { + (n / v).ln() - 1.0 + } + })); eos.lift().residual_helmholtz_energy(t, v, &n) / t + ig }, D::from(0.0), diff --git a/crates/feos/tests/pcsaft/critical_point.rs b/crates/feos/tests/pcsaft/critical_point.rs index 6e4ee3156..9183e5e5d 100644 --- a/crates/feos/tests/pcsaft/critical_point.rs +++ b/crates/feos/tests/pcsaft/critical_point.rs @@ -1,10 +1,11 @@ use approx::assert_relative_eq; use feos::pcsaft::{PcSaft, PcSaftParameters}; -use feos_core::State; use feos_core::parameter::IdentifierOption; +use feos_core::{SolverOptions, State}; use nalgebra::dvector; use quantity::*; use std::error::Error; +use std::sync::Arc; #[test] fn test_critical_point_pure() -> Result<(), Box> { @@ -46,3 +47,31 @@ fn test_critical_point_mix() -> Result<(), Box> { ); Ok(()) } + +#[test] +fn test_critical_point_limits() -> Result<(), Box> { + let params = PcSaftParameters::from_json( + vec!["propane", "butane"], + "tests/pcsaft/test_parameters.json", + None, + IdentifierOption::Name, + )?; + let options = SolverOptions { + verbosity: feos_core::Verbosity::Iter, + ..Default::default() + }; + let saft = Arc::new(PcSaft::new(params)); + let cp_pure = State::critical_point_pure(&saft, None, None, options)?; + println!("{} {}", cp_pure[0], cp_pure[1]); + let molefracs = dvector![0.0, 1.0]; + let cp_2 = State::critical_point(&saft, Some(&molefracs), None, None, options)?; + println!("{}", cp_2); + let molefracs = dvector![1.0, 0.0]; + let cp_1 = State::critical_point(&saft, Some(&molefracs), None, None, options)?; + println!("{}", cp_1); + assert_eq!(cp_pure[0].temperature, cp_1.temperature); + assert_eq!(cp_pure[0].density, cp_1.density); + assert_eq!(cp_pure[1].temperature, cp_2.temperature); + assert_eq!(cp_pure[1].density, cp_2.density); + Ok(()) +} diff --git a/crates/feos/tests/pcsaft/tp_flash.rs b/crates/feos/tests/pcsaft/tp_flash.rs index ffe45e0b4..b31a14998 100644 --- a/crates/feos/tests/pcsaft/tp_flash.rs +++ b/crates/feos/tests/pcsaft/tp_flash.rs @@ -4,7 +4,6 @@ use feos_core::parameter::IdentifierOption; use feos_core::{Contributions, FeosResult, PhaseEquilibrium, SolverOptions}; use nalgebra::dvector; use quantity::*; -use std::error::Error; fn read_params(components: Vec<&str>) -> FeosResult { PcSaftParameters::from_json( @@ -16,7 +15,7 @@ fn read_params(components: Vec<&str>) -> FeosResult { } #[test] -fn test_tp_flash() -> Result<(), Box> { +fn test_tp_flash() -> FeosResult<()> { let propane = PcSaft::new(read_params(vec!["propane"])?); let butane = PcSaft::new(read_params(vec!["butane"])?); let t = 250.0 * KELVIN; @@ -63,3 +62,51 @@ fn test_tp_flash() -> Result<(), Box> { ); Ok(()) } + +#[test] +fn test_tp_flash_zero_component() -> FeosResult<()> { + let eos_full = PcSaft::new(PcSaftParameters::from_json( + vec!["propane", "butane", "hexane"], + "tests/pcsaft/test_parameters.json", + None, + IdentifierOption::Name, + )?); + let eos_binary = PcSaft::new(PcSaftParameters::from_json( + vec!["butane", "hexane"], + "tests/pcsaft/test_parameters.json", + None, + IdentifierOption::Name, + )?); + let options = SolverOptions { + verbosity: feos_core::Verbosity::Iter, + ..Default::default() + }; + let vle_full = PhaseEquilibrium::tp_flash( + &&eos_full, + 300.0 * KELVIN, + 1.2 * BAR, + &(dvector![0.0, 0.5, 0.5] * MOL), + None, + options, + None, + )?; + let vle_binary = PhaseEquilibrium::tp_flash( + &&eos_binary, + 300.0 * KELVIN, + 1.2 * BAR, + &(dvector![0.5, 0.5] * MOL), + None, + options, + None, + )?; + println!("{vle_full}\n{vle_binary}"); + assert_eq!( + vle_full.liquid().molefracs[1], + vle_binary.liquid().molefracs[0] + ); + assert_eq!( + vle_full.vapor().molefracs[1], + vle_binary.vapor().molefracs[0] + ); + Ok(()) +}