From 7bc4e43119375828fcad4a20d66787ed1e70e3ec Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Fri, 17 Apr 2026 16:33:22 -0400 Subject: [PATCH] Update gas mixture with correct n calculation --- src/mixtures/gas_mixture.rs | 85 ++++++++++++++++++++++++++++++++-- src/properties/mod.rs | 2 + src/properties/test_helpers.rs | 45 ++++++++++++++++++ src/properties/thermo_fit.rs | 5 ++ 4 files changed, 132 insertions(+), 5 deletions(-) create mode 100644 src/properties/test_helpers.rs diff --git a/src/mixtures/gas_mixture.rs b/src/mixtures/gas_mixture.rs index f0df99f..f302a01 100644 --- a/src/mixtures/gas_mixture.rs +++ b/src/mixtures/gas_mixture.rs @@ -1,19 +1,62 @@ +use std::collections::HashMap; + use crate::{ consts::P_REF, - properties::{ - thermo_fit::{Phase, SpeciesThermoData}, - transport_fit::SpeciesTransportData, - }, + matrix::Matrix, + properties::thermo_fit::{Phase, SpeciesThermoData}, }; pub struct GasMixture { pub(crate) ns: Vec, pub(crate) nsum: f64, pub(crate) species: Vec, - pub(crate) transport_data: Vec, + pub(crate) coeffs: Matrix, + pub(crate) elements: HashMap, } impl GasMixture { + pub fn new(ns_kg_moles: &[f64], species: &[SpeciesThermoData]) -> Self { + // First calculate the total and per species kg-moles/kg mixture + let (nsum_kg_moles, mass_sum) = ns_kg_moles + .iter() + .zip(species.iter()) + .fold((0.0, 0.0), |acc, (n, s)| { + (acc.0 + n, acc.1 + n * s.molecular_weight) + }); + + let nsum = nsum_kg_moles / mass_sum; + let ns = ns_kg_moles.iter().map(|n| n / mass_sum).collect(); + + // Now build out the element list + let mut elements = HashMap::new(); + let mut ei = 0; + for s in species.iter() { + for element in s.elements.iter() { + if !elements.contains_key(&element.element) { + elements.insert(element.element.clone(), ei); + ei += 1 + } + } + } + + // 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); + }); + } + + GasMixture { + ns, + nsum, + species: species.to_vec(), + coeffs, + elements, + } + } // 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 { @@ -78,3 +121,35 @@ impl GasMixture { .collect() } } + +#[cfg(test)] +mod test { + use crate::{ + assert_delta, assert_vec_delta, mixtures::gas_mixture::GasMixture, + properties::test_helpers::h2_o2_h2o_thermo_data, + }; + + #[test] + pub fn test_gas_mixture_creation() { + let species = h2_o2_h2o_thermo_data(); + + let ns = [1.0, 2.0, 3.0]; + + let gas = GasMixture::new(&ns, &species); + + let expected_ns = [ + 0.008329211761713246, + 0.01665842352342649, + 0.02498763528513974, + ]; + assert_vec_delta!(expected_ns, gas.ns, 1e-12); + assert_delta!(gas.nsum, 0.04997527057027948, 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); + } +} diff --git a/src/properties/mod.rs b/src/properties/mod.rs index f160f70..c177906 100644 --- a/src/properties/mod.rs +++ b/src/properties/mod.rs @@ -5,5 +5,7 @@ pub mod thermo_fit; pub mod transport_db; pub mod transport_fit; mod utils; +#[cfg(test)] +pub mod test_helpers; pub use error::PropertiesError; diff --git a/src/properties/test_helpers.rs b/src/properties/test_helpers.rs new file mode 100644 index 0000000..3f6d093 --- /dev/null +++ b/src/properties/test_helpers.rs @@ -0,0 +1,45 @@ +use super::thermo_fit::{Phase, SpeciesElement, SpeciesThermoData, ThermoPolynomial}; + +pub fn simple_thermo_polynomial( + cp_over_r: f64, + h_b: f64, + s_b: f64, + temp_range: (f64, f64), +) -> ThermoPolynomial { + ThermoPolynomial { + a: vec![0.0, 0.0, cp_over_r, 0.0, 0.0, 0.0, 0.0, h_b, s_b], + temp_range, + } +} + +pub fn h2_o2_h2o_thermo_data() -> Vec { + vec![ + SpeciesThermoData::new( + "H2", + vec![SpeciesElement { element: "H".to_string(), count: 2.0 }], + Phase::Gas, + vec![simple_thermo_polynomial(3.5, -1012.0, 29.11, (200.0, 6000.0))], + 2.01594, + 0.0, + ), + SpeciesThermoData::new( + "O2", + vec![SpeciesElement { element: "O".to_string(), count: 2.0 }], + Phase::Gas, + vec![simple_thermo_polynomial(3.5, -1005.0, 30.03, (200.0, 6000.0))], + 31.9988, + 0.0, + ), + SpeciesThermoData::new( + "H2O", + vec![ + SpeciesElement { element: "H".to_string(), count: 2.0 }, + SpeciesElement { element: "O".to_string(), count: 1.0 }, + ], + Phase::Gas, + vec![simple_thermo_polynomial(4.0, -29192.0, 23.03, (200.0, 6000.0))], + 18.01528, + -241826.0, + ), + ] +} diff --git a/src/properties/thermo_fit.rs b/src/properties/thermo_fit.rs index 893bbcf..272263a 100644 --- a/src/properties/thermo_fit.rs +++ b/src/properties/thermo_fit.rs @@ -1,8 +1,10 @@ +#[derive(Clone)] pub enum Phase { Gas, Condensed, } +#[derive(Clone)] pub struct SpeciesThermoData { pub name: String, pub elements: Vec, @@ -11,11 +13,14 @@ pub struct SpeciesThermoData { pub molecular_weight: f64, pub h_formation: f64, } + +#[derive(Clone)] pub struct ThermoPolynomial { pub a: Vec, pub temp_range: (f64, f64), } +#[derive(Clone)] pub struct SpeciesElement { pub element: String, pub count: f64,