use super::parameters::{PcSaftAssociationRecord, PcSaftParameters, PcSaftPars};
use crate::association::{Association, AssociationStrength};
use crate::hard_sphere::{HardSphere, HardSphereProperties, MonomerShape};
use feos_core::{
EntropyScaling, Molarweight, ReferenceSystem, Residual, ResidualDyn, StateHD, Subset,
};
use nalgebra::DVector;
use num_dual::{DualNum, partial2};
use quantity::ad::first_derivative;
use quantity::*;
use std::f64::consts::{FRAC_PI_6, PI};
use typenum::P2;
pub(crate) mod dispersion;
pub(crate) mod hard_chain;
mod pcsaft_binary;
mod pcsaft_pure;
pub(crate) mod polar;
use dispersion::Dispersion;
use hard_chain::HardChain;
pub use pcsaft_binary::PcSaftBinary;
pub use pcsaft_pure::PcSaftPure;
pub use polar::DQVariants;
use polar::{Dipole, DipoleQuadrupole, Quadrupole};
#[derive(Copy, Clone)]
pub struct PcSaftOptions {
pub max_eta: f64,
pub max_iter_cross_assoc: usize,
pub tol_cross_assoc: f64,
pub dq_variant: DQVariants,
}
impl Default for PcSaftOptions {
fn default() -> Self {
Self {
max_eta: 0.5,
max_iter_cross_assoc: 50,
tol_cross_assoc: 1e-10,
dq_variant: DQVariants::DQ35,
}
}
}
pub struct PcSaft {
pub parameters: PcSaftParameters,
params: PcSaftPars,
options: PcSaftOptions,
hard_chain: bool,
dipole: bool,
quadrupole: bool,
dipole_quadrupole: bool,
association: Option<Association>,
}
impl PcSaft {
pub fn new(parameters: PcSaftParameters) -> Self {
Self::with_options(parameters, PcSaftOptions::default())
}
pub fn with_options(parameters: PcSaftParameters, options: PcSaftOptions) -> Self {
let params = PcSaftPars::new(¶meters);
let hard_chain = params.m.iter().any(|m| (m - 1.0).abs() > 1e-15);
let dipole = params.ndipole > 0;
let quadrupole = params.nquadpole > 0;
let dipole_quadrupole = params.ndipole > 0 && params.nquadpole > 0;
let association = (!parameters.association.is_empty())
.then(|| Association::new(options.max_iter_cross_assoc, options.tol_cross_assoc));
Self {
parameters,
params,
options,
hard_chain,
dipole,
quadrupole,
dipole_quadrupole,
association,
}
}
}
impl ResidualDyn for PcSaft {
fn components(&self) -> usize {
self.parameters.pure.len()
}
fn compute_max_density<D: DualNum<f64> + Copy>(&self, molefracs: &DVector<D>) -> D {
let msigma3 = self
.params
.m
.component_mul(&self.params.sigma.map(|v| v.powi(3)));
(msigma3.map(D::from).dot(molefracs) * FRAC_PI_6).recip() * self.options.max_eta
}
fn reduced_helmholtz_energy_density_contributions<D: DualNum<f64> + Copy>(
&self,
state: &StateHD<D>,
) -> Vec<(&'static str, D)> {
let mut v = Vec::with_capacity(7);
let d = self.params.hs_diameter(state.temperature);
v.push((
"Hard Sphere",
HardSphere.helmholtz_energy_density(&self.params, state),
));
if self.hard_chain {
v.push((
"Hard Chain",
HardChain.helmholtz_energy_density(&self.params, state),
))
}
v.push((
"Dispersion",
Dispersion.helmholtz_energy_density(&self.params, state),
));
if self.dipole {
v.push((
"Dipole",
Dipole.helmholtz_energy_density(&self.params, state),
))
}
if self.quadrupole {
v.push((
"Quadrupole",
Quadrupole.helmholtz_energy_density(&self.params, state),
))
}
if self.dipole_quadrupole {
v.push((
"DipoleQuadrupole",
DipoleQuadrupole.helmholtz_energy_density(
&self.params,
state,
self.options.dq_variant,
),
))
}
if let Some(association) = self.association.as_ref() {
v.push((
"Association",
association.helmholtz_energy_density(
&self.params,
&self.parameters.association,
state,
&d,
),
))
}
v
}
}
impl Subset for PcSaft {
fn subset(&self, component_list: &[usize]) -> Self {
Self::with_options(self.parameters.subset(component_list), self.options)
}
}
impl Molarweight for PcSaft {
fn molar_weight(&self) -> MolarWeight<DVector<f64>> {
self.parameters.molar_weight.clone()
}
}
impl HardSphereProperties for PcSaftPars {
fn monomer_shape<N: DualNum<f64>>(&self, _: N) -> MonomerShape<'_, N> {
MonomerShape::NonSpherical(self.m.map(N::from))
}
fn hs_diameter<D: DualNum<f64> + Copy>(&self, temperature: D) -> DVector<D> {
let ti = temperature.recip() * -3.0;
DVector::from_fn(self.sigma.len(), |i, _| {
-((ti * self.epsilon_k[i]).exp() * 0.12 - 1.0) * self.sigma[i]
})
}
}
impl AssociationStrength for PcSaftPars {
type Record = PcSaftAssociationRecord;
fn association_strength_ij<D: DualNum<f64> + Copy>(
&self,
temperature: D,
comp_i: usize,
comp_j: usize,
association_parameters_ij: &Self::Record,
) -> D {
let f_ab = (temperature.recip() * association_parameters_ij.epsilon_k_ab).exp_m1();
let k_ab = association_parameters_ij.kappa_ab
* (self.sigma[comp_i] * self.sigma[comp_j]).powf(1.5);
f_ab * k_ab
}
}
fn omega11(t: f64) -> f64 {
1.06036 * t.powf(-0.15610)
+ 0.19300 * (-0.47635 * t).exp()
+ 1.03587 * (-1.52996 * t).exp()
+ 1.76474 * (-3.89411 * t).exp()
}
fn omega22(t: f64) -> f64 {
1.16145 * t.powf(-0.14874) + 0.52487 * (-0.77320 * t).exp() + 2.16178 * (-2.43787 * t).exp()
- 6.435e-4 * t.powf(0.14874) * (18.0323 * t.powf(-0.76830) - 7.27371).sin()
}
#[inline]
fn chapman_enskog_thermal_conductivity(
temperature: Temperature,
molarweight: MolarWeight,
m: f64,
sigma: f64,
epsilon_k: f64,
) -> ThermalConductivity {
let t = temperature.into_reduced();
0.083235 * (t * m / molarweight.convert_to(GRAM / MOL)).sqrt()
/ sigma.powi(2)
/ omega22(t / epsilon_k)
* WATT
/ METER
/ KELVIN
}
impl EntropyScaling for PcSaft {
fn viscosity_reference(
&self,
temperature: Temperature,
_: Volume,
moles: &Moles<DVector<f64>>,
) -> Viscosity {
let p = &self.params;
let mw = &self.parameters.molar_weight;
let x = (moles / moles.sum()).into_value();
let ce: Vec<_> = (0..self.components())
.map(|i| {
let tr = (temperature / p.epsilon_k[i] / KELVIN).into_value();
5.0 / 16.0 * (mw.get(i) * KB / NAV * temperature / PI).sqrt()
/ omega22(tr)
/ (p.sigma[i] * ANGSTROM).powi::<P2>()
})
.collect();
let mut ce_mix = 0.0 * MILLI * PASCAL * SECOND;
for i in 0..self.components() {
let denom: f64 = (0..self.components())
.map(|j| {
x[j] * (1.0
+ (ce[i] / ce[j]).into_value().sqrt()
* (mw.get(j) / mw.get(i)).powf(1.0 / 4.0))
.powi(2)
/ (8.0 * (1.0 + (mw.get(i) / mw.get(j)).into_value())).sqrt()
})
.sum();
ce_mix += ce[i] * x[i] / denom
}
ce_mix
}
fn viscosity_correlation(&self, s_res: f64, x: &DVector<f64>) -> f64 {
let coefficients = self
.params
.viscosity
.as_ref()
.expect("Missing viscosity coefficients.");
let m = x.dot(&self.params.m);
let s = s_res / m;
let pref = x.component_mul(&self.params.m) / m;
let a = coefficients.row(0).transpose().dot(x);
let b = coefficients.row(1).transpose().dot(&pref);
let c = coefficients.row(2).transpose().dot(&pref);
let d = coefficients.row(3).transpose().dot(&pref);
a + b * s + c * s.powi(2) + d * s.powi(3)
}
fn diffusion_reference(
&self,
temperature: Temperature,
volume: Volume,
moles: &Moles<DVector<f64>>,
) -> Diffusivity {
if self.components() != 1 {
panic!("Diffusion coefficients in PC-SAFT are only implemented for pure components!");
}
let p = &self.params;
let mw = &self.parameters.molar_weight;
let density = moles.sum() / volume;
let res: Vec<_> = (0..self.components())
.map(|i| {
let tr = (temperature / p.epsilon_k[i] / KELVIN).into_value();
3.0 / 8.0 / (p.sigma[i] * ANGSTROM).powi::<P2>() / omega11(tr) / (density * NAV)
* (temperature * RGAS / PI / mw.get(i) / p.m[i]).sqrt()
})
.collect();
res[0]
}
fn diffusion_correlation(&self, s_res: f64, x: &DVector<f64>) -> f64 {
if self.components() != 1 {
panic!("Diffusion coefficients in PC-SAFT are only implemented for pure components!");
}
let coefficients = self
.params
.diffusion
.as_ref()
.expect("Missing diffusion coefficients.");
let m = x.dot(&self.params.m);
let s = s_res / m;
let pref = x.component_mul(&self.params.m) / m;
let a = coefficients.row(0).transpose().dot(x);
let b = coefficients.row(1).transpose().dot(&pref);
let c = coefficients.row(2).transpose().dot(&pref);
let d = coefficients.row(3).transpose().dot(&pref);
let e = coefficients.row(4).transpose().dot(&pref);
a + b * s - c * (1.0 - s.exp()) * s.powi(2) - d * s.powi(4) - e * s.powi(8)
}
fn thermal_conductivity_reference(
&self,
temperature: Temperature,
volume: Volume,
moles: &Moles<DVector<f64>>,
) -> ThermalConductivity {
if self.components() != 1 {
panic!("Thermal conductivity in PC-SAFT is only implemented for pure components!");
}
let p = &self.params;
let mws = self.molar_weight();
let (_, s_res) = first_derivative(
partial2(
|t, &v, n| -self.residual_helmholtz_energy_unit(t, v, n) / n.sum(),
&volume,
moles,
),
temperature,
);
let res: Vec<_> = (0..self.components())
.map(|i| {
let tr = (temperature / p.epsilon_k[i] / KELVIN).into_value();
let s_res_reduced = s_res.into_reduced() / p.m[i];
let ref_ce = chapman_enskog_thermal_conductivity(
temperature,
mws.get(i),
p.m[i],
p.sigma[i],
p.epsilon_k[i],
);
let alpha_visc = (-s_res_reduced / -0.5).exp();
let ref_ts = (-0.0167141 * tr / p.m[i] + 0.0470581 * (tr / p.m[i]).powi(2))
* (p.m[i] * p.m[i] * p.sigma[i].powi(3) * p.epsilon_k[i])
* 1e-5
* WATT
/ METER
/ KELVIN;
ref_ce + ref_ts * alpha_visc
})
.collect();
res[0]
}
fn thermal_conductivity_correlation(&self, s_res: f64, x: &DVector<f64>) -> f64 {
if self.components() != 1 {
panic!("Thermal conductivity in PC-SAFT is only implemented for pure components!");
}
let coefficients = self
.params
.thermal_conductivity
.as_ref()
.expect("Missing thermal conductivity coefficients");
let m = x.dot(&self.params.m);
let s = s_res / m;
let pref = x.component_mul(&self.params.m) / m;
let a = coefficients.row(0).transpose().dot(x);
let b = coefficients.row(1).transpose().dot(&pref);
let c = coefficients.row(2).transpose().dot(&pref);
let d = coefficients.row(3).transpose().dot(&pref);
a + b * s + c * (1.0 - s.exp()) + d * s.powi(2)
}
}
#[cfg(test)]
mod tests {
use super::*;
use crate::pcsaft::parameters::utils::{
butane_parameters, nonane_heptane_parameters, nonane_parameters, propane_butane_parameters,
propane_parameters,
};
use approx::assert_relative_eq;
use feos_core::*;
use nalgebra::dvector;
use quantity::{BAR, KELVIN, METER, PASCAL, RGAS};
use typenum::{P2, P3};
#[test]
fn ideal_gas_pressure() {
let e = &propane_parameters();
let t = 200.0 * KELVIN;
let v = 1e-3 * METER.powi::<P3>();
let n = dvector![1.0] * MOL;
let s = State::new_nvt(&e, t, v, &n).unwrap();
let p_ig = s.total_moles * RGAS * t / v;
assert_relative_eq!(s.pressure(Contributions::IdealGas), p_ig, epsilon = 1e-10);
assert_relative_eq!(
s.pressure(Contributions::IdealGas) + s.pressure(Contributions::Residual),
s.pressure(Contributions::Total),
epsilon = 1e-10
);
}
#[test]
fn ideal_gas_heat_capacity_joback() {
let e = &propane_parameters();
let t = 200.0 * KELVIN;
let v = 1e-3 * METER.powi::<P3>();
let n = dvector![1.0] * MOL;
let s = State::new_nvt(&e, t, v, &n).unwrap();
let p_ig = s.total_moles * RGAS * t / v;
assert_relative_eq!(s.pressure(Contributions::IdealGas), p_ig, epsilon = 1e-10);
assert_relative_eq!(
s.pressure(Contributions::IdealGas) + s.pressure(Contributions::Residual),
s.pressure(Contributions::Total),
epsilon = 1e-10
);
}
#[test]
fn hard_sphere() {
let t = 250.0;
let v = 1000.0;
let n = 1.0;
let s = StateHD::new(t, v, &dvector![n]);
let a_rust = HardSphere
.helmholtz_energy_density(&propane_parameters().params as &PcSaftPars, &s)
* v;
assert_relative_eq!(a_rust, 0.410610492598808, epsilon = 1e-10);
}
#[test]
fn hard_sphere_mix() {
let t = 250.0;
let v = 1000.0;
let n = 1.0;
let s = StateHD::new(t, v, &dvector![n]);
let a1 =
HardSphere.helmholtz_energy_density(&propane_parameters().params as &PcSaftPars, &s);
let a2 =
HardSphere.helmholtz_energy_density(&butane_parameters().params as &PcSaftPars, &s);
let s1m = StateHD::new(t, v, &dvector![n, 0.0]);
let a1m = HardSphere
.helmholtz_energy_density(&propane_butane_parameters().params as &PcSaftPars, &s1m);
let s2m = StateHD::new(t, v, &dvector![0.0, n]);
let a2m = HardSphere
.helmholtz_energy_density(&propane_butane_parameters().params as &PcSaftPars, &s2m);
assert_relative_eq!(a1, a1m, epsilon = 1e-14);
assert_relative_eq!(a2, a2m, epsilon = 1e-14);
}
#[test]
fn new_tpn() {
let e = &propane_parameters();
let t = 300.0 * KELVIN;
let p = BAR;
let m = dvector![1.5] * MOL;
let s = State::new_npt(&e, t, p, &m, None);
let p_calc = if let Ok(state) = s {
state.pressure(Contributions::Total)
} else {
0.0 * PASCAL
};
assert_relative_eq!(p, p_calc, epsilon = 1e-6);
}
#[test]
fn vle_pure() {
let e = &propane_parameters();
let t = 300.0 * KELVIN;
let vle = PhaseEquilibrium::pure(&e, t, None, Default::default());
if let Ok(v) = vle {
assert_relative_eq!(
v.vapor().pressure(Contributions::Total),
v.liquid().pressure(Contributions::Total),
epsilon = 1e-6
)
}
}
#[test]
fn critical_point() {
let e = &propane_parameters();
let t = 300.0 * KELVIN;
let cp = State::critical_point(&e, None, Some(t), None, Default::default());
if let Ok(v) = cp {
assert_relative_eq!(v.temperature, 375.1244078318015 * KELVIN, epsilon = 1e-8)
}
}
#[test]
fn mix_single() {
let e1 = &propane_parameters();
let e2 = &butane_parameters();
let e12 = &propane_butane_parameters();
let t = 300.0 * KELVIN;
let v = 0.02456883872966545 * METER.powi::<P3>();
let m1 = dvector![2.0] * MOL;
let m1m = dvector![2.0, 0.0] * MOL;
let m2m = dvector![0.0, 2.0] * MOL;
let s1 = State::new_nvt(&e1, t, v, &m1).unwrap();
let s2 = State::new_nvt(&e2, t, v, &m1).unwrap();
let s1m = State::new_nvt(&e12, t, v, &m1m).unwrap();
let s2m = State::new_nvt(&e12, t, v, &m2m).unwrap();
assert_relative_eq!(
s1.pressure(Contributions::Total),
s1m.pressure(Contributions::Total),
epsilon = 1e-12
);
assert_relative_eq!(
s2.pressure(Contributions::Total),
s2m.pressure(Contributions::Total),
epsilon = 1e-12
)
}
#[test]
fn viscosity() -> FeosResult<()> {
let e = &propane_parameters();
let t = 300.0 * KELVIN;
let p = BAR;
let n = dvector![1.0] * MOL;
let s = State::new_npt(&e, t, p, &n, None)?;
assert_relative_eq!(
s.viscosity(),
0.00797 * MILLI * PASCAL * SECOND,
epsilon = 1e-5
);
assert_relative_eq!(
s.ln_viscosity_reduced(),
(s.viscosity() / e.viscosity_reference(s.temperature, s.volume, &s.moles))
.into_value()
.ln(),
epsilon = 1e-15
);
Ok(())
}
#[test]
fn viscosity_mix() -> FeosResult<()> {
let e = &nonane_heptane_parameters();
let nonane = &nonane_parameters();
let t = 303.15 * KELVIN;
let p = 500.0 * BAR;
let n = dvector![0.25, 0.75] * MOL;
let viscosity_mix = State::new_npt(&e, t, p, &n, None)?.viscosity();
let viscosity_paper = 0.68298 * MILLI * PASCAL * SECOND;
assert_relative_eq!(viscosity_paper, viscosity_mix, epsilon = 1e-8);
let n_pseudo_mix = dvector![1.0, 0.0] * MOL;
let viscosity_pseudo_mix = State::new_npt(&e, t, p, &n_pseudo_mix, None)?.viscosity();
let n_nonane = dvector![1.0] * MOL;
let viscosity_nonane = State::new_npt(&nonane, t, p, &n_nonane, None)?.viscosity();
assert_relative_eq!(viscosity_pseudo_mix, viscosity_nonane, epsilon = 1e-15);
Ok(())
}
#[test]
fn diffusion() -> FeosResult<()> {
let e = &propane_parameters();
let t = 300.0 * KELVIN;
let p = BAR;
let n = dvector![1.0] * MOL;
let s = State::new_npt(&e, t, p, &n, None)?;
assert_relative_eq!(
s.diffusion(),
0.01505 * (CENTI * METER).powi::<P2>() / SECOND,
epsilon = 1e-5
);
assert_relative_eq!(
s.ln_diffusion_reduced(),
(s.diffusion() / e.diffusion_reference(s.temperature, s.volume, &s.moles))
.into_value()
.ln(),
epsilon = 1e-15
);
Ok(())
}
}
#[cfg(test)]
mod tests_parameter_fit {
use super::pcsaft_binary::test::pcsaft_binary;
use super::pcsaft_pure::test::pcsaft;
use super::*;
use approx::assert_relative_eq;
use feos_core::DensityInitialization::Liquid;
use feos_core::{Contributions, PropertiesAD, ReferenceSystem};
use feos_core::{FeosResult, ParametersAD, PhaseEquilibrium, State};
use nalgebra::{U1, U3, U8, vector};
use num_dual::{DualStruct, DualVec, partial};
use quantity::{BAR, KELVIN, LITER, MOL, PASCAL};
fn pcsaft_non_assoc() -> PcSaftPure<f64, 4> {
let m = 1.5;
let sigma = 3.4;
let epsilon_k = 180.0;
let mu = 2.2;
let params = [m, sigma, epsilon_k, mu];
PcSaftPure(params)
}
#[test]
fn test_vapor_pressure_derivatives() -> FeosResult<()> {
let pcsaft_params = [
"m",
"sigma",
"epsilon_k",
"mu",
"kappa_ab",
"epsilon_k_ab",
"na",
"nb",
];
let (pcsaft, _) = pcsaft()?;
let pcsaft_ad = pcsaft.named_derivatives(pcsaft_params);
let temperature = 250.0 * KELVIN;
let p = pcsaft_ad.vapor_pressure(temperature)?;
let p = p.convert_into(PASCAL);
let (p, grad) = (p.re, p.eps.unwrap_generic(U8, U1));
println!("{p:.5}");
println!("{grad:.5?}");
for (i, par) in pcsaft_params.into_iter().enumerate() {
let mut params = pcsaft.0;
let h = params[i] * 1e-7;
params[i] += h;
let pcsaft_h = PcSaftPure(params);
let (p_h, _) =
PhaseEquilibrium::pure_t(&pcsaft_h, temperature, None, Default::default())?;
let dp_h = (p_h.convert_into(PASCAL) - p) / h;
let dp = grad[i];
println!(
"{par:12}: {:11.5} {:11.5} {:.3e}",
dp_h,
dp,
((dp_h - dp) / dp).abs()
);
assert_relative_eq!(dp, dp_h, max_relative = 1e-6);
}
Ok(())
}
#[test]
fn test_vapor_pressure_derivatives_fit() -> FeosResult<()> {
let pcsaft = pcsaft_non_assoc();
let pcsaft_ad = pcsaft.named_derivatives(["m", "sigma", "epsilon_k"]);
let temperature = 150.0 * KELVIN;
let p = pcsaft_ad.vapor_pressure(temperature)?;
let p = p.convert_into(PASCAL);
let (p, grad) = (p.re, p.eps.unwrap_generic(U3, U1));
println!("{p:.5}");
println!("{grad:.5?}");
for (i, par) in ["m", "sigma", "epsilon_k"].into_iter().enumerate() {
let mut params = pcsaft.0;
let h = params[i] * 1e-7;
params[i] += h;
let pcsaft_h = PcSaftPure(params);
let (p_h, _) =
PhaseEquilibrium::pure_t(&pcsaft_h, temperature, None, Default::default())?;
let dp_h = (p_h.convert_into(PASCAL) - p) / h;
let dp = grad[i];
println!(
"{par:12}: {:11.5} {:11.5} {:.3e}",
dp_h,
dp,
((dp_h - dp) / dp).abs()
);
assert_relative_eq!(dp, dp_h, max_relative = 1e-6);
}
Ok(())
}
#[test]
fn test_equilibrium_liquid_density_derivatives_fit() -> FeosResult<()> {
let pcsaft = pcsaft_non_assoc();
let pcsaft_ad = pcsaft.named_derivatives(["m", "sigma", "epsilon_k"]);
let temperature = 150.0 * KELVIN;
let (p, rho) = pcsaft_ad.equilibrium_liquid_density(temperature)?;
let p = p.convert_into(PASCAL);
let rho = rho.convert_into(MOL / LITER);
let (p, p_grad) = (p.re, p.eps.unwrap_generic(U3, U1));
let (rho, rho_grad) = (rho.re, rho.eps.unwrap_generic(U3, U1));
println!("{p:.5} {rho:.5}");
println!("{p_grad:.5?}");
println!("{rho_grad:.5?}");
for (i, par) in ["m", "sigma", "epsilon_k"].into_iter().enumerate() {
let mut params = pcsaft.0;
let h = params[i] * 1e-7;
params[i] += h;
let pcsaft_h = PcSaftPure(params);
let (p_h, [_, rho_h]) =
PhaseEquilibrium::pure_t(&pcsaft_h, temperature, None, Default::default())?;
let dp_h = (p_h.convert_into(PASCAL) - p) / h;
let drho_h = (rho_h.convert_into(MOL / LITER) - rho) / h;
let dp = p_grad[i];
let drho = rho_grad[i];
println!(
"{par:12}: {:11.5} {:11.5} {:.3e} {:11.5} {:11.5} {:.3e}",
dp_h,
dp,
((dp_h - dp) / dp).abs(),
drho_h,
drho,
((drho_h - drho) / drho).abs()
);
assert_relative_eq!(dp, dp_h, max_relative = 1e-6);
}
Ok(())
}
#[test]
fn test_liquid_density_derivatives_fit() -> FeosResult<()> {
let pcsaft = pcsaft_non_assoc();
let pcsaft_ad = pcsaft.named_derivatives(["m", "sigma", "epsilon_k"]);
let temperature = 150.0 * KELVIN;
let pressure = BAR;
let rho = pcsaft_ad.liquid_density(temperature, pressure)?;
let rho = rho.convert_into(MOL / LITER);
let (rho, grad) = (rho.re, rho.eps.unwrap_generic(U3, U1));
println!("{rho:.5}");
println!("{grad:.5?}");
for (i, par) in ["m", "sigma", "epsilon_k"].into_iter().enumerate() {
let mut params = pcsaft.0;
let h = params[i] * 1e-7;
params[i] += h;
let pcsaft_h = PcSaftPure(params);
let rho_h = State::new_xpt(
&pcsaft_h,
temperature,
pressure,
&vector![1.0],
Some(Liquid),
)?
.density;
let drho_h = (rho_h.convert_into(MOL / LITER) - rho) / h;
let drho = grad[i];
println!(
"{par:12}: {:11.5} {:11.5} {:.3e}",
drho_h,
drho,
((drho_h - drho) / drho).abs()
);
assert_relative_eq!(drho, drho_h, max_relative = 1e-6);
}
Ok(())
}
#[test]
fn test_bubble_point_pressure() -> FeosResult<()> {
let (pcsaft, _) = pcsaft_binary()?;
let pcsaft_ad = pcsaft.named_derivatives(["k_ij"]);
let temperature = 500.0 * KELVIN;
let x = vector![0.5, 0.5];
let p = pcsaft_ad.bubble_point_pressure(temperature, None, x)?;
let p = p.convert_into(BAR);
let (p, [[grad]]) = (p.re, p.eps.unwrap_generic(U1, U1).data.0);
println!("{p:.5}");
println!("{grad:.5?}");
let (params, mut kij) = pcsaft.0;
let h = 1e-7;
kij += h;
let pcsaft_h = PcSaftBinary::new(params, kij);
let p_h = PhaseEquilibrium::bubble_point(
&pcsaft_h,
temperature,
&x,
None,
None,
Default::default(),
)?
.vapor()
.pressure(Contributions::Total);
let dp_h = (p_h.convert_into(BAR) - p) / h;
println!(
"k_ij: {:11.5} {:11.5} {:.3e}",
dp_h,
grad,
((dp_h - grad) / grad).abs()
);
assert_relative_eq!(grad, dp_h, max_relative = 1e-6);
Ok(())
}
#[test]
fn test_dew_point_pressure() -> FeosResult<()> {
let (pcsaft, _) = pcsaft_binary()?;
let pcsaft_ad = pcsaft.named_derivatives(["k_ij"]);
let temperature = 500.0 * KELVIN;
let y = vector![0.5, 0.5];
let p = pcsaft_ad.dew_point_pressure(temperature, None, y)?;
let p = p.convert_into(BAR);
let (p, [[grad]]) = (p.re, p.eps.unwrap_generic(U1, U1).data.0);
println!("{p:.5}");
println!("{grad:.5?}");
let (params, mut kij) = pcsaft.0;
let h = 1e-7;
kij += h;
let pcsaft_h = PcSaftBinary::new(params, kij);
let p_h = PhaseEquilibrium::dew_point(
&pcsaft_h,
temperature,
&y,
None,
None,
Default::default(),
)?
.vapor()
.pressure(Contributions::Total);
let dp_h = (p_h.convert_into(BAR) - p) / h;
println!(
"k_ij: {:11.5} {:11.5} {:.3e}",
dp_h,
grad,
((dp_h - grad) / grad).abs()
);
assert_relative_eq!(grad, dp_h, max_relative = 1e-6);
Ok(())
}
#[test]
fn test_bubble_point_temperature() -> FeosResult<()> {
let (pcsaft, _) = pcsaft_binary()?;
let pcsaft_ad = pcsaft.named_derivatives(["k_ij"]);
let pressure = Pressure::from_reduced(DualVec::from(45. * BAR.into_reduced()));
let t_init = Temperature::from_reduced(DualVec::from(500.0));
let x = vector![0.5, 0.5].map(DualVec::from);
let t = PhaseEquilibrium::bubble_point(
&pcsaft_ad,
pressure,
&x,
Some(t_init),
None,
Default::default(),
)?
.vapor()
.temperature;
let t = t.convert_into(KELVIN);
let (t, [[grad]]) = (t.re, t.eps.unwrap_generic(U1, U1).data.0);
println!("{t:.5}");
println!("{grad:.5?}");
let (params, mut kij) = pcsaft.0;
let h = 1e-7;
kij += h;
let pcsaft_h = PcSaftBinary::new(params, kij);
let t_h = PhaseEquilibrium::bubble_point(
&pcsaft_h,
pressure.re(),
&x.map(|x| x.re()),
Some(t_init.re()),
None,
Default::default(),
)?
.vapor()
.temperature;
let dt_h = (t_h.convert_into(KELVIN) - t) / h;
println!(
"k_ij: {:11.5} {:11.5} {:.3e}",
dt_h,
grad,
((dt_h - grad) / grad).abs()
);
assert_relative_eq!(grad, dt_h, max_relative = 1e-6);
Ok(())
}
#[test]
fn test_dew_point_temperature() -> FeosResult<()> {
let (pcsaft, _) = pcsaft_binary()?;
let pcsaft_ad = pcsaft.named_derivatives(["k_ij"]);
let pressure = Pressure::from_reduced(DualVec::from(45. * BAR.into_reduced()));
let t_init = Temperature::from_reduced(DualVec::from(500.0));
let x = vector![0.5, 0.5].map(DualVec::from);
let t = PhaseEquilibrium::dew_point(
&pcsaft_ad,
pressure,
&x,
Some(t_init),
None,
Default::default(),
)?
.vapor()
.temperature;
let t = t.convert_into(KELVIN);
let (t, [[grad]]) = (t.re, t.eps.unwrap_generic(U1, U1).data.0);
println!("{t:.5}");
println!("{grad:.5?}");
let (params, mut kij) = pcsaft.0;
let h = 1e-7;
kij += h;
let pcsaft_h = PcSaftBinary::new(params, kij);
let t_h = PhaseEquilibrium::dew_point(
&pcsaft_h,
pressure.re(),
&x.map(|x| x.re()),
Some(t_init.re()),
None,
Default::default(),
)?
.vapor()
.temperature;
let dt_h = (t_h.convert_into(KELVIN) - t) / h;
println!(
"k_ij: {:11.5} {:11.5} {:.3e}",
dt_h,
grad,
((dt_h - grad) / grad).abs()
);
assert_relative_eq!(grad, dt_h, max_relative = 1e-6);
Ok(())
}
#[test]
fn test_bubble_point_temperature_derivative() -> FeosResult<()> {
let (pcsaft, _) = pcsaft_binary()?;
let pcsaft_ad = pcsaft;
let mut temperature = 500.0 * KELVIN;
let x = vector![0.5, 0.5];
let (p, grad) = first_derivative(
partial(
|t, x| {
let eos = pcsaft_ad.lift();
PhaseEquilibrium::bubble_point(&eos, t, x, None, None, Default::default())
.unwrap()
.vapor()
.pressure(Contributions::Total)
},
&x,
),
temperature,
);
println!("{p:.5}");
println!("{grad:.5?}");
let h = 1e-6 * KELVIN;
temperature += h;
let p_h = PhaseEquilibrium::bubble_point(
&pcsaft_ad,
temperature,
&x,
None,
None,
Default::default(),
)?
.vapor()
.pressure(Contributions::Total);
let dp_h = (p_h - p) / h;
println!(
"{:11.5?} {:11.5?} {:.3e}",
dp_h,
grad,
((dp_h - grad).convert_into(grad)).abs()
);
assert_relative_eq!(grad, dp_h, max_relative = 1e-7);
Ok(())
}
#[test]
fn test_bubble_point_pressure_derivative() -> FeosResult<()> {
let (pcsaft, _) = pcsaft_binary()?;
let pcsaft_ad = pcsaft;
let mut pressure = 45. * BAR;
let t0 = Some(500. * KELVIN);
let x = vector![0.5, 0.5];
let (t, grad) = first_derivative(
partial2(
|p, x, &t0| {
let eos = pcsaft_ad.lift();
PhaseEquilibrium::bubble_point(&eos, p, x, t0, None, Default::default())
.unwrap()
.vapor()
.temperature
},
&x,
&t0,
),
pressure,
);
println!("{t:.5}");
println!("{grad:.5?}");
let h = 1e-5 * BAR;
pressure += h;
let t_h =
PhaseEquilibrium::bubble_point(&pcsaft_ad, pressure, &x, t0, None, Default::default())?
.vapor()
.temperature;
let dt_h = (t_h - t) / h;
println!(
"{:11.5?} {:11.5?} {:.3e}",
dt_h,
grad,
((dt_h - grad).convert_into(grad)).abs()
);
assert_relative_eq!(grad, dt_h, max_relative = 1e-7);
Ok(())
}
}