diff --git a/src/mixtures/gas_mixture.rs b/src/mixtures/gas_mixture.rs index 65085eb..d88eeb4 100644 --- a/src/mixtures/gas_mixture.rs +++ b/src/mixtures/gas_mixture.rs @@ -6,18 +6,11 @@ use crate::{ properties::thermo_fit::{Phase, SpeciesThermoData}, }; -pub struct MixtureComponent { - //kg-moles component/kg_mixture - n: f64, - // Coefficients - a: Vec, - s: SpeciesThermoData, -} - pub struct GasMixture { + pub(crate) ns: Vec, pub(crate) nsum: f64, - pub(crate) gasses: Vec, - pub(crate) condensed: Vec, + pub(crate) species: Vec, + pub(crate) coeffs: Matrix, pub(crate) elements: HashMap, pub(crate) binitial: Vec, } @@ -47,30 +40,23 @@ impl GasMixture { } } + // Now build the coefficients + let mut coeffs = Matrix::new(ei, species.len(), 0.0); + for (j, s) in species.iter().enumerate() { + s.elements.iter().for_each(|e| { + // Safe to unwrap because elements should always have e + let i = *elements.get(&e.element).unwrap(); + coeffs.set(i, j, e.count); + }); + } + let binitial = get_b_current(&elements, species, &ns); - // Now separate SpeciesThermoData in gas and condensed MixtureComponents - let (gasses, condensed) = ns - .iter() - .zip(species.iter()) - .map(|(n, s)| { - let mut a = vec![0.0; elements.len()]; - s.elements.iter().for_each(|e| { - let i = elements.get(&e.element).unwrap(); - a[*i] += e.count; - }); - MixtureComponent { - n: *n, - a, - s: s.clone(), - } - }) - .partition(|c| matches!(c.s.phase, Phase::Gas)); - GasMixture { + ns, nsum, - gasses, - condensed, + species: species.to_vec(), + coeffs, elements, binitial, } @@ -78,18 +64,18 @@ impl GasMixture { // 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 + self.ns .iter() - .chain(self.condensed.iter()) - .map(|c| -> f64 { - match c.s.phase { + .zip(self.species.iter()) + .map(|(n, s)| -> f64 { + match s.phase { Phase::Gas => { - let p = - c.s.polynomial_at(temp) - .expect("Gas doesn't have a polynomial"); + let p = 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() + + (n / self.nsum).ln() } Phase::Condensed => todo!(), } @@ -101,16 +87,16 @@ impl GasMixture { // // Equations 2.17 from reference paper pub fn gas_entropies_over_rt(&self, temp: f64, pressure: f64) -> Vec { - self.gasses + self.ns .iter() - .chain(self.condensed.iter()) - .map(|c| -> f64 { - match c.s.phase { + .zip(self.species.iter()) + .map(|(n, s)| -> f64 { + match 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() + let p = s + .polynomial_at(temp) + .expect("Gas doesn't have a polynomial"); + p.s_over_r(temp) - (n / self.nsum).ln() - (pressure / P_REF).ln() } Phase::Condensed => todo!(), } @@ -122,16 +108,16 @@ impl GasMixture { // 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 + self.ns .iter() - .chain(self.condensed.iter()) - .map(|c| -> f64 { - match c.s.phase { + .zip(self.species.iter()) + .map(|(n, s)| -> f64 { + match 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) + let p = s + .polynomial_at(temp) + .expect("Gas doesn't have a polynomial"); + n * p.h_over_rt(temp) } Phase::Condensed => todo!(), } @@ -177,17 +163,15 @@ mod test { 0.01665842352342649, 0.02498763528513974, ]; - let ns: Vec = gas.gasses.iter().map(|c| c.n).collect(); - - assert_vec_delta!(expected_ns, ns, 1e-12); + assert_vec_delta!(expected_ns, gas.ns, 1e-12); assert_delta!(gas.nsum, 0.04997527057027948, 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); + assert_delta!(gas.coeffs.get(0, 0).unwrap(), 2.0, 1e-12); + assert_delta!(gas.coeffs.get(0, 1).unwrap(), 0.0, 1e-12); + assert_delta!(gas.coeffs.get(0, 2).unwrap(), 2.0, 1e-12); + assert_delta!(gas.coeffs.get(1, 0).unwrap(), 0.0, 1e-12); + assert_delta!(gas.coeffs.get(1, 1).unwrap(), 2.0, 1e-12); + assert_delta!(gas.coeffs.get(1, 2).unwrap(), 1.0, 1e-12); let expected_b = [0.06663369409370597, 0.05830448233199272]; assert_vec_delta!(gas.binitial, expected_b, 1e-12);