Skip to content
This repository was archived by the owner on Jun 14, 2022. It is now read-only.

Pr impl residual #51

Merged
prehner merged 6 commits into from
Apr 6, 2022
Merged

Pr impl residual #51

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
07fc8b2
Add loss function to estimator, and Antoine equation for vapor pressu…
g-bauer Feb 25, 2022
23a21ff
Fixed due to renaming
g-bauer Mar 14, 2022
2860d16
resolve merge conflict
g-bauer Mar 22, 2022
7a1bad4
Merge branch 'main' of github.com:feos-org/feos-core
g-bauer Apr 5, 2022
4993953
Impl PR Helmholtz energy as contribution
g-bauer Mar 22, 2022
17c0fd3
updated docstrings
g-bauer Apr 5, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
82 changes: 46 additions & 36 deletions src/cubic.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@ use num_dual::DualNum;
use quantity::si::{SIArray1, SIUnit};
use serde::{Deserialize, Serialize};
use std::f64::consts::SQRT_2;
use std::fmt;
use std::rc::Rc;

const KB_A3: f64 = 13806490.0;
Expand Down Expand Up @@ -165,12 +166,51 @@ impl Parameter for PengRobinsonParameters {
}
}

struct PengRobinsonContribution {
parameters: Rc<PengRobinsonParameters>,
}

impl<D: DualNum<f64>> HelmholtzEnergyDual<D> for PengRobinsonContribution {
fn helmholtz_energy(&self, state: &StateHD<D>) -> D {
// temperature dependent a parameter
let p = &self.parameters;
let x = &state.molefracs;
let ak = (&p.tc.mapv(|tc| (D::one() - (state.temperature / tc).sqrt())) * &p.kappa + 1.0)
.mapv(|x| x.powi(2))
* &p.a;

// Mixing rules
let mut ak_mix = D::zero();
for i in 0..ak.len() {
for j in 0..ak.len() {
ak_mix += (ak[i] * ak[j]).sqrt() * (x[i] * x[j] * (1.0 - p.k_ij[(i, j)]));
}
}
let b = (x * &p.b).sum();

// Helmholtz energy
let n = state.moles.sum();
let v = state.volume;
n * ((v / (v - b * n)).ln()
- ak_mix / (b * SQRT_2 * 2.0 * state.temperature)
* ((v * (SQRT_2 - 1.0) + b * n) / (v * (SQRT_2 + 1.0) - b * n)).ln())
}
}

impl fmt::Display for PengRobinsonContribution {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
write!(f, "Peng Robinson")
}
}

/// A simple version of the Peng-Robinson equation of state.
pub struct PengRobinson {
/// Parameters
parameters: Rc<PengRobinsonParameters>,
/// Ideal gas contributions to the Helmholtz energy
ideal_gas: Joback,
/// Non-ideal contributions to the Helmholtz energy
contributions: Vec<Box<dyn HelmholtzEnergy>>,
}

impl PengRobinson {
Expand All @@ -180,9 +220,14 @@ impl PengRobinson {
|| Joback::default(parameters.tc.len()),
|j| Joback::new(j.clone()),
);
let contributions: Vec<Box<dyn HelmholtzEnergy>> =
vec![Box::new(PengRobinsonContribution {
parameters: parameters.clone(),
})];
Self {
parameters,
ideal_gas,
contributions,
}
}
}
Expand All @@ -202,42 +247,7 @@ impl EquationOfState for PengRobinson {
}

fn residual(&self) -> &[Box<dyn HelmholtzEnergy>] {
unreachable!()
}

fn evaluate_residual<D: DualNum<f64>>(&self, state: &StateHD<D>) -> D {
// temperature dependent a parameter
let p = &self.parameters;
let x = &state.molefracs;
let ak = (&p.tc.mapv(|tc| (D::one() - (state.temperature / tc).sqrt())) * &p.kappa + 1.0)
.mapv(|x| x.powi(2))
* &p.a;

// Mixing rules
let mut ak_mix = D::zero();
for i in 0..ak.len() {
for j in 0..ak.len() {
ak_mix += (ak[i] * ak[j]).sqrt() * (x[i] * x[j] * (1.0 - p.k_ij[(i, j)]));
}
}
let b = (x * &p.b).sum();

// Helmholtz energy
let n = state.moles.sum();
let v = state.volume;
n * ((v / (v - b * n)).ln()
- ak_mix / (b * SQRT_2 * 2.0 * state.temperature)
* ((v * (SQRT_2 - 1.0) + b * n) / (v * (SQRT_2 + 1.0) - b * n)).ln())
}

fn evaluate_residual_contributions<D: DualNum<f64>>(
&self,
state: &StateHD<D>,
) -> Vec<(String, D)>
where
dyn HelmholtzEnergy: HelmholtzEnergyDual<D>,
{
vec![("Peng-Robinson".into(), self.evaluate_residual(state))]
&self.contributions
}

fn ideal_gas(&self) -> &dyn IdealGasContribution {
Expand Down
8 changes: 1 addition & 7 deletions src/equation_of_state.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -165,9 +165,6 @@ pub trait EquationOfState {
fn residual(&self) -> &[Box<dyn HelmholtzEnergy>];

/// Evaluate the residual reduced Helmholtz energy $\beta A^\mathrm{res}$.
///
/// For simple equations of state (see e.g. `PengRobinson`) it might be
/// easier to overwrite this function instead of implementing `residual`.
fn evaluate_residual<D: DualNum<f64>>(&self, state: &StateHD<D>) -> D
where
dyn HelmholtzEnergy: HelmholtzEnergyDual<D>,
Expand All @@ -179,10 +176,7 @@ pub trait EquationOfState {
}

/// Evaluate the reduced Helmholtz energy of each individual contribution
/// and return them together with a string representatino of the contribution.
///
/// If `evaluate_residual` is implemented instead of `residual`, this function
/// also needs to be overwritten to avoid panics.
/// and return them together with a string representation of the contribution.
fn evaluate_residual_contributions<D: DualNum<f64>>(
&self,
state: &StateHD<D>,
Expand Down