diff --git a/src/mixtures/gas_mixture.rs b/src/mixtures/gas_mixture.rs index 3d98367..65085eb 100644 --- a/src/mixtures/gas_mixture.rs +++ b/src/mixtures/gas_mixture.rs @@ -2,15 +2,16 @@ use std::collections::HashMap; use crate::{ consts::P_REF, + matrix::Matrix, properties::thermo_fit::{Phase, SpeciesThermoData}, }; pub struct MixtureComponent { //kg-moles component/kg_mixture - pub(crate) n: f64, + n: f64, // Coefficients - pub(crate) coeff: Vec, - pub(crate) s: SpeciesThermoData, + a: Vec, + s: SpeciesThermoData, } pub struct GasMixture { @@ -21,39 +22,6 @@ pub struct GasMixture { pub(crate) binitial: Vec, } -impl MixtureComponent { - // µ/RT for the current mixture component at fixed temp and pressure - // µ = Molar Gibbs free energy = H - TS - // Therefore µ/RT = H/RT - S/R - pub fn chem_potential_over_rt(&self, temp: f64, pressure: f64, nsum: f64) -> f64 { - match self.s.phase { - Phase::Gas => { - let p = self - .s - .polynomial_at(temp) - .expect("Gas doesn't have a polynomial"); - p.h_over_rt(temp) - p.s_over_r(temp) - + (pressure / P_REF).ln() - + (self.n / nsum).ln() - } - Phase::Condensed => todo!(), - } - } - - pub fn entropy_over_r(&self, temp: f64, pressure: f64, nsum: f64) -> f64 { - match self.s.phase { - Phase::Gas => { - let p = self - .s - .polynomial_at(temp) - .expect("Gas doesn't have a polynomial"); - p.s_over_r(temp) - (self.n / nsum).ln() - (pressure / P_REF).ln() - } - Phase::Condensed => todo!(), - } - } -} - impl GasMixture { pub fn new(ns_kg_moles: &[f64], species: &[SpeciesThermoData]) -> Self { // First calculate the total and per species kg-moles/kg mixture @@ -79,8 +47,10 @@ impl GasMixture { } } + let binitial = get_b_current(&elements, species, &ns); + // Now separate SpeciesThermoData in gas and condensed MixtureComponents - let (gasses, condensed): (Vec<_>, Vec<_>) = ns + let (gasses, condensed) = ns .iter() .zip(species.iter()) .map(|(n, s)| { @@ -91,14 +61,12 @@ impl GasMixture { }); MixtureComponent { n: *n, - coeff: a, + a, s: s.clone(), } }) .partition(|c| matches!(c.s.phase, Phase::Gas)); - let binitial = get_b_current(&elements, gasses.iter().chain(condensed.iter())); - GasMixture { nsum, gasses, @@ -107,26 +75,83 @@ impl GasMixture { binitial, } } + // Calculate the normalized chemical potential (μ/RT) for each component in the mixture. + // Equations 2.11 from reference paper + pub fn gas_chem_potentials_over_rt(&self, temp: f64, pressure: f64) -> Vec { + self.gasses + .iter() + .chain(self.condensed.iter()) + .map(|c| -> f64 { + match c.s.phase { + Phase::Gas => { + let p = + c.s.polynomial_at(temp) + .expect("Gas doesn't have a polynomial"); + p.h_over_rt(temp) - p.s_over_r(temp) + + (pressure / P_REF).ln() + + (c.n / self.nsum).ln() + } + Phase::Condensed => todo!(), + } + }) + .collect() + } - pub fn get_b_current(&self) -> Vec { - get_b_current( - &self.elements, - self.gasses.iter().chain(self.condensed.iter()), - ) + // Calculate the normalized entropy (S/R) for each mixture component + // + // Equations 2.17 from reference paper + pub fn gas_entropies_over_rt(&self, temp: f64, pressure: f64) -> Vec { + self.gasses + .iter() + .chain(self.condensed.iter()) + .map(|c| -> f64 { + match c.s.phase { + Phase::Gas => { + let p = + c.s.polynomial_at(temp) + .expect("Gas doesn't have a polynomial"); + p.s_over_r(temp) - (c.n / self.nsum).ln() - (pressure / P_REF).ln() + } + Phase::Condensed => todo!(), + } + }) + .collect() + } + + // Calculate the normalized mixture enthalpy (H/RT) + // Note that the enthalpy doesn't have a dependence on the pressure. + // Equation 2.14 from the paper + pub fn mixture_h_over_rt(&self, temp: f64) -> Vec { + self.gasses + .iter() + .chain(self.condensed.iter()) + .map(|c| -> f64 { + match c.s.phase { + Phase::Gas => { + let p = + c.s.polynomial_at(temp) + .expect("Gas doesn't have a polynomial"); + c.n * p.h_over_rt(temp) + } + Phase::Condensed => todo!(), + } + }) + .collect() } } // Current kilogram-atoms of element i per kg mixture -pub fn get_b_current<'a>( +pub fn get_b_current( ele_map: &HashMap, - components: impl IntoIterator, + species: &[SpeciesThermoData], + ns: &[f64], ) -> Vec { let mut b = vec![0.0; ele_map.len()]; - for c in components { - for e in c.s.elements.iter() { + for (n, s) in ns.iter().zip(species) { + for e in s.elements.iter() { let i = *ele_map.get(&e.element).unwrap(); - b[i] += c.n * e.count; + b[i] += n * e.count; } } b @@ -157,12 +182,12 @@ mod test { assert_vec_delta!(expected_ns, ns, 1e-12); assert_delta!(gas.nsum, 0.04997527057027948, 1e-12); - assert_delta!(gas.gasses[0].coeff[0], 2.0, 1e-12); - assert_delta!(gas.gasses[1].coeff[0], 0.0, 1e-12); - assert_delta!(gas.gasses[2].coeff[0], 2.0, 1e-12); - assert_delta!(gas.gasses[0].coeff[1], 0.0, 1e-12); - assert_delta!(gas.gasses[1].coeff[1], 2.0, 1e-12); - assert_delta!(gas.gasses[2].coeff[1], 1.0, 1e-12); + assert_delta!(gas.gasses[0].a[0], 2.0, 1e-12); + assert_delta!(gas.gasses[1].a[0], 0.0, 1e-12); + assert_delta!(gas.gasses[2].a[0], 2.0, 1e-12); + assert_delta!(gas.gasses[0].a[1], 0.0, 1e-12); + assert_delta!(gas.gasses[1].a[1], 2.0, 1e-12); + assert_delta!(gas.gasses[2].a[1], 1.0, 1e-12); let expected_b = [0.06663369409370597, 0.05830448233199272]; assert_vec_delta!(gas.binitial, expected_b, 1e-12); diff --git a/src/solvers/equations.rs b/src/solvers/equations.rs index bf50720..0829845 100644 --- a/src/solvers/equations.rs +++ b/src/solvers/equations.rs @@ -35,11 +35,11 @@ impl std::error::Error for SolverError {} // Solve Ax=b for x using Gauss-Jordan elimination // A must be a square n x n matrix and b must be a column vector, i.e. n x 1 // Algorithm taken from Numerical Methods for Engineers by Ayyub and McCuen -pub fn gauss_jordan_elimination(a: &Matrix, y: &Matrix) -> Result, CEAError> { - validate_ax_eq_b(a, y)?; +pub fn gauss_jordan_elimination(a: &Matrix, b: &Matrix) -> Result, CEAError> { + validate_ax_eq_b(a, b)?; let mut a = a.clone(); - let mut b = y.clone(); + let mut b = b.clone(); for i in 0..a.shape().0 { // First pivot rows