mirror of
https://github.com/aselimov/cea-rs.git
synced 2026-04-19 00:24:20 +00:00
GitHub
commit
44e2426aa2
4 changed files with 444 additions and 31 deletions
|
|
@ -98,14 +98,14 @@ fn parse_species<'a>(
|
|||
lines.next().ok_or(PropertiesError::InvalidFile)?;
|
||||
}
|
||||
|
||||
Ok(SpeciesThermoData {
|
||||
name,
|
||||
polynomials,
|
||||
Ok(SpeciesThermoData::new(
|
||||
&name,
|
||||
elements,
|
||||
phase,
|
||||
polynomials,
|
||||
molecular_weight,
|
||||
h_formation,
|
||||
})
|
||||
))
|
||||
}
|
||||
|
||||
fn parse_polynomials_block<'a>(
|
||||
|
|
@ -295,9 +295,17 @@ mod test {
|
|||
-1.742171366e+01,
|
||||
];
|
||||
|
||||
assert_vec_delta!(species.polynomials[0].a, real_coeff_1, 1e-9);
|
||||
assert_delta!(species.polynomials[0].temp_range.0, 300.000, 1e-3);
|
||||
assert_delta!(species.polynomials[0].temp_range.1, 1000.000, 1e-3);
|
||||
assert_vec_delta!(species.polynomial_at(650.0).unwrap().a, real_coeff_1, 1e-9);
|
||||
assert_delta!(
|
||||
species.polynomial_at(650.0).unwrap().temp_range.0,
|
||||
300.000,
|
||||
1e-3
|
||||
);
|
||||
assert_delta!(
|
||||
species.polynomial_at(650.0).unwrap().temp_range.1,
|
||||
1000.000,
|
||||
1e-3
|
||||
);
|
||||
|
||||
let real_coeff_2 = [
|
||||
-3.523782900e+05,
|
||||
|
|
@ -311,9 +319,17 @@ mod test {
|
|||
-2.695610360e+00,
|
||||
];
|
||||
|
||||
assert_vec_delta!(species.polynomials[1].a, real_coeff_2, 1e-9);
|
||||
assert_delta!(species.polynomials[1].temp_range.0, 1000.000, 1e-3);
|
||||
assert_delta!(species.polynomials[1].temp_range.1, 6000.000, 1e-3);
|
||||
assert_vec_delta!(species.polynomial_at(3500.0).unwrap().a, real_coeff_2, 1e-9);
|
||||
assert_delta!(
|
||||
species.polynomial_at(3500.0).unwrap().temp_range.0,
|
||||
1000.000,
|
||||
1e-3
|
||||
);
|
||||
assert_delta!(
|
||||
species.polynomial_at(3500.0).unwrap().temp_range.1,
|
||||
6000.000,
|
||||
1e-3
|
||||
);
|
||||
}
|
||||
|
||||
#[test]
|
||||
|
|
@ -363,10 +379,10 @@ END REACTANTS
|
|||
assert!(matches!(alcl3.phase, Phase::Gas));
|
||||
assert_delta!(alcl3.molecular_weight, 133.3405380, 1e-7);
|
||||
assert_delta!(alcl3.h_formation, -584678.863, 1e-3);
|
||||
assert_eq!(alcl3.polynomials.len(), 2);
|
||||
assert_eq!(alcl3.num_polynomials(), 2);
|
||||
|
||||
assert_vec_delta!(
|
||||
alcl3.polynomials[0].a,
|
||||
alcl3.polynomial_at(650.0).unwrap().a,
|
||||
[
|
||||
7.750600970e+04,
|
||||
-1.440779717e+03,
|
||||
|
|
@ -380,11 +396,19 @@ END REACTANTS
|
|||
],
|
||||
1e-9
|
||||
);
|
||||
assert_delta!(alcl3.polynomials[0].temp_range.0, 300.0, 1e-3);
|
||||
assert_delta!(alcl3.polynomials[0].temp_range.1, 1000.0, 1e-3);
|
||||
assert_delta!(
|
||||
alcl3.polynomial_at(650.0).unwrap().temp_range.0,
|
||||
300.0,
|
||||
1e-3
|
||||
);
|
||||
assert_delta!(
|
||||
alcl3.polynomial_at(650.0).unwrap().temp_range.1,
|
||||
1000.0,
|
||||
1e-3
|
||||
);
|
||||
|
||||
assert_vec_delta!(
|
||||
alcl3.polynomials[1].a,
|
||||
alcl3.polynomial_at(3500.0).unwrap().a,
|
||||
[
|
||||
-1.378630916e+05,
|
||||
-5.579207290e+01,
|
||||
|
|
@ -398,8 +422,16 @@ END REACTANTS
|
|||
],
|
||||
1e-9
|
||||
);
|
||||
assert_delta!(alcl3.polynomials[1].temp_range.0, 1000.0, 1e-3);
|
||||
assert_delta!(alcl3.polynomials[1].temp_range.1, 6000.0, 1e-3);
|
||||
assert_delta!(
|
||||
alcl3.polynomial_at(3500.0).unwrap().temp_range.0,
|
||||
1000.0,
|
||||
1e-3
|
||||
);
|
||||
assert_delta!(
|
||||
alcl3.polynomial_at(3500.0).unwrap().temp_range.1,
|
||||
6000.0,
|
||||
1e-3
|
||||
);
|
||||
|
||||
// --- Air (reactant 0) ---
|
||||
let air = &thermo_db.reactants[0];
|
||||
|
|
@ -416,10 +448,10 @@ END REACTANTS
|
|||
assert!(matches!(air.phase, Phase::Gas));
|
||||
assert_delta!(air.molecular_weight, 28.9651159, 1e-7);
|
||||
assert_delta!(air.h_formation, -125.530, 1e-3);
|
||||
assert_eq!(air.polynomials.len(), 2);
|
||||
assert_eq!(air.num_polynomials(), 2);
|
||||
|
||||
assert_vec_delta!(
|
||||
air.polynomials[0].a,
|
||||
air.polynomial_at(650.0).unwrap().a,
|
||||
[
|
||||
1.009950160e+04,
|
||||
-1.968275610e+02,
|
||||
|
|
@ -433,11 +465,11 @@ END REACTANTS
|
|||
],
|
||||
1e-9
|
||||
);
|
||||
assert_delta!(air.polynomials[0].temp_range.0, 300.0, 1e-3);
|
||||
assert_delta!(air.polynomials[0].temp_range.1, 1000.0, 1e-3);
|
||||
assert_delta!(air.polynomial_at(650.0).unwrap().temp_range.0, 300.0, 1e-3);
|
||||
assert_delta!(air.polynomial_at(650.0).unwrap().temp_range.1, 1000.0, 1e-3);
|
||||
|
||||
assert_vec_delta!(
|
||||
air.polynomials[1].a,
|
||||
air.polynomial_at(3500.0).unwrap().a,
|
||||
[
|
||||
2.415214430e+05,
|
||||
-1.257874600e+03,
|
||||
|
|
@ -451,8 +483,16 @@ END REACTANTS
|
|||
],
|
||||
1e-9
|
||||
);
|
||||
assert_delta!(air.polynomials[1].temp_range.0, 1000.0, 1e-3);
|
||||
assert_delta!(air.polynomials[1].temp_range.1, 6000.0, 1e-3);
|
||||
assert_delta!(
|
||||
air.polynomial_at(3500.0).unwrap().temp_range.0,
|
||||
1000.0,
|
||||
1e-3
|
||||
);
|
||||
assert_delta!(
|
||||
air.polynomial_at(3500.0).unwrap().temp_range.1,
|
||||
6000.0,
|
||||
1e-3
|
||||
);
|
||||
|
||||
// --- n-Butanol (reactant 1) ---
|
||||
let butanol = &thermo_db.reactants[1];
|
||||
|
|
@ -467,6 +507,6 @@ END REACTANTS
|
|||
assert!(matches!(butanol.phase, Phase::Condensed));
|
||||
assert_delta!(butanol.molecular_weight, 74.1216000, 1e-7);
|
||||
assert_delta!(butanol.h_formation, -278510.000, 1e-3);
|
||||
assert_eq!(butanol.polynomials.len(), 0);
|
||||
assert_eq!(butanol.num_polynomials(), 0);
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -7,7 +7,7 @@ pub struct SpeciesThermoData {
|
|||
pub name: String,
|
||||
pub elements: Vec<SpeciesElement>,
|
||||
pub phase: Phase,
|
||||
pub polynomials: Vec<ThermoPolynomial>,
|
||||
polynomials: Vec<ThermoPolynomial>,
|
||||
pub molecular_weight: f64,
|
||||
pub h_formation: f64,
|
||||
}
|
||||
|
|
@ -20,3 +20,193 @@ pub struct SpeciesElement {
|
|||
pub element: String,
|
||||
pub count: f64,
|
||||
}
|
||||
|
||||
impl SpeciesThermoData {
|
||||
pub fn new(
|
||||
name: &str,
|
||||
elements: Vec<SpeciesElement>,
|
||||
phase: Phase,
|
||||
polynomials: Vec<ThermoPolynomial>,
|
||||
molecular_weight: f64,
|
||||
h_formation: f64,
|
||||
) -> Self {
|
||||
Self {
|
||||
name: name.to_string(),
|
||||
elements,
|
||||
phase,
|
||||
polynomials,
|
||||
molecular_weight,
|
||||
h_formation,
|
||||
}
|
||||
}
|
||||
|
||||
pub fn num_polynomials(&self) -> usize {
|
||||
self.polynomials.len()
|
||||
}
|
||||
|
||||
pub fn polynomial_at(&self, temp: f64) -> Option<&ThermoPolynomial> {
|
||||
//TODO: Not the most efficient. Can refactor to pre-compute tables
|
||||
//and do 1-d linear interpolation if needed
|
||||
//
|
||||
//TODO: I Think condensed species need to be treated differently. Verify how that works in
|
||||
//the paper.
|
||||
if self.polynomials.is_empty() {
|
||||
return None;
|
||||
}
|
||||
|
||||
let i_polynomial = self
|
||||
.polynomials
|
||||
.iter()
|
||||
.rposition(|polynomial| temp > polynomial.temp_range.0)
|
||||
.unwrap_or(0);
|
||||
Some(&self.polynomials[i_polynomial])
|
||||
}
|
||||
}
|
||||
|
||||
impl ThermoPolynomial {
|
||||
/// Calculate using eq 4.9 from reference paper
|
||||
/// NOTE: This is normalized and unitless
|
||||
pub fn cp_over_r(&self, temp: f64) -> f64 {
|
||||
let inv_temp = 1.0 / temp;
|
||||
self.a[0] * inv_temp * inv_temp
|
||||
+ self.a[1] * inv_temp
|
||||
+ self.a[2]
|
||||
+ self.a[3] * temp
|
||||
+ self.a[4] * temp * temp
|
||||
+ self.a[5] * temp * temp * temp
|
||||
+ self.a[6] * temp * temp * temp * temp
|
||||
}
|
||||
/// Calculate using eq 4.10 from reference paper
|
||||
/// NOTE: This is normalized and unitless
|
||||
pub fn h_over_rt(&self, temp: f64) -> f64 {
|
||||
let inv_temp = 1.0 / temp;
|
||||
-self.a[0] * inv_temp * inv_temp
|
||||
+ self.a[1] * inv_temp * temp.ln()
|
||||
+ self.a[2]
|
||||
+ self.a[3] * temp / 2.0
|
||||
+ self.a[4] * temp * temp / 3.0
|
||||
+ self.a[5] * temp * temp * temp / 4.0
|
||||
+ self.a[6] * temp * temp * temp * temp / 5.0
|
||||
+ self.a[7] * inv_temp
|
||||
}
|
||||
/// Calculate using eq 4.11 from reference paper
|
||||
/// NOTE: This is normalized and unitless
|
||||
pub fn s_over_r(&self, temp: f64) -> f64 {
|
||||
let inv_temp = 1.0 / temp;
|
||||
-self.a[0] * inv_temp * inv_temp / 2.0 - self.a[1] * inv_temp
|
||||
+ self.a[2] * temp.ln()
|
||||
+ self.a[3] * temp
|
||||
+ self.a[4] * temp * temp / 2.0
|
||||
+ self.a[5] * temp * temp * temp / 3.0
|
||||
+ self.a[6] * temp * temp * temp * temp / 4.0
|
||||
+ self.a[8]
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod test {
|
||||
use super::ThermoPolynomial;
|
||||
use crate::{
|
||||
assert_delta,
|
||||
properties::thermo_fit::{Phase, SpeciesThermoData},
|
||||
};
|
||||
|
||||
#[test]
|
||||
fn test_cp_over_r() {
|
||||
let result = ThermoPolynomial {
|
||||
a: vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0],
|
||||
temp_range: (0.0, 0.0),
|
||||
}
|
||||
.cp_over_r(100.0);
|
||||
assert_delta!(result, 706050403.0201, 1e-4);
|
||||
|
||||
let result = ThermoPolynomial {
|
||||
a: vec![4.0, 2.0, 1.0, 2.0, 1.0, 1.0, 1.0],
|
||||
temp_range: (0.0, 0.0),
|
||||
}
|
||||
.cp_over_r(2.0);
|
||||
assert_delta!(result, 35.0, 1e-10);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_h_over_rt() {
|
||||
let result = ThermoPolynomial {
|
||||
a: vec![1.0, 2.0, 3.0, 12.0, 9.0, 8.0, 5.0, 8.0],
|
||||
temp_range: (0.0, 0.0),
|
||||
}
|
||||
.h_over_rt(2.0);
|
||||
assert_delta!(result, 63.44314718055995, 1e-10);
|
||||
|
||||
let result = ThermoPolynomial {
|
||||
a: vec![4.0, 0.0, 3.0, 4.0, 0.0, 0.0, 0.0, 2.0],
|
||||
temp_range: (0.0, 0.0),
|
||||
}
|
||||
.h_over_rt(100.0);
|
||||
assert_delta!(result, 203.0196, 1e-4);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_s_over_r() {
|
||||
let result = ThermoPolynomial {
|
||||
a: vec![2.0, 3.0, 4.0, 1.0, 5.0, 2.0, 8.0, 1.0, 12.0],
|
||||
temp_range: (0.0, 0.0),
|
||||
}
|
||||
.s_over_r(100.0);
|
||||
assert_delta!(result, 200691797.0572474, 1e-7);
|
||||
|
||||
let result = ThermoPolynomial {
|
||||
a: vec![4.0, 2.0, 0.0, 2.0, 2.0, 0.0, 0.0, 0.0, 3.0],
|
||||
temp_range: (0.0, 0.0),
|
||||
}
|
||||
.s_over_r(2.0);
|
||||
assert_delta!(result, 9.5, 1e-10);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_polynomial_at() {
|
||||
let polynomials = vec![
|
||||
ThermoPolynomial {
|
||||
a: vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0],
|
||||
temp_range: (1.0, 100.0),
|
||||
},
|
||||
ThermoPolynomial {
|
||||
a: vec![4.0, 2.0, 1.0, 2.0, 1.0, 1.0, 1.0],
|
||||
temp_range: (100.0, 300.0),
|
||||
},
|
||||
];
|
||||
|
||||
let data = SpeciesThermoData {
|
||||
name: "".to_string(),
|
||||
elements: vec![],
|
||||
phase: Phase::Gas,
|
||||
polynomials,
|
||||
molecular_weight: 1.0,
|
||||
h_formation: 1.0,
|
||||
};
|
||||
|
||||
assert!(std::ptr::eq(
|
||||
data.polynomial_at(0.5).unwrap(),
|
||||
&data.polynomials[0]
|
||||
));
|
||||
assert!(std::ptr::eq(
|
||||
data.polynomial_at(50.0).unwrap(),
|
||||
&data.polynomials[0]
|
||||
));
|
||||
assert!(std::ptr::eq(
|
||||
data.polynomial_at(100.0).unwrap(),
|
||||
&data.polynomials[0]
|
||||
));
|
||||
assert!(std::ptr::eq(
|
||||
data.polynomial_at(100.0 + 1e-12).unwrap(),
|
||||
&data.polynomials[1]
|
||||
));
|
||||
assert!(std::ptr::eq(
|
||||
data.polynomial_at(150.0 + 1e-12).unwrap(),
|
||||
&data.polynomials[1]
|
||||
));
|
||||
assert!(std::ptr::eq(
|
||||
data.polynomial_at(500.0 + 1e-12).unwrap(),
|
||||
&data.polynomials[1]
|
||||
));
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -50,11 +50,11 @@ fn parse_species_transport_block<'a>(
|
|||
(ViscosityOrConductivity::Conductivity, fit) => conductivities.push(fit),
|
||||
}
|
||||
}
|
||||
Ok(SpeciesTransportData {
|
||||
name,
|
||||
Ok(SpeciesTransportData::new(
|
||||
&name,
|
||||
viscosities,
|
||||
conductivities,
|
||||
})
|
||||
))
|
||||
}
|
||||
|
||||
fn parse_species_header_line(line: &str) -> Result<(String, usize, usize), PropertiesError> {
|
||||
|
|
|
|||
|
|
@ -1,7 +1,45 @@
|
|||
pub struct SpeciesTransportData {
|
||||
pub name: String,
|
||||
pub viscosities: Vec<TransportFit>,
|
||||
pub conductivities: Vec<TransportFit>,
|
||||
pub(crate) viscosities: Vec<TransportFit>,
|
||||
pub(crate) conductivities: Vec<TransportFit>,
|
||||
}
|
||||
|
||||
impl SpeciesTransportData {
|
||||
pub fn new(
|
||||
name: &str,
|
||||
viscosities: Vec<TransportFit>,
|
||||
conductivities: Vec<TransportFit>,
|
||||
) -> Self {
|
||||
SpeciesTransportData {
|
||||
name: name.to_string(),
|
||||
viscosities,
|
||||
conductivities,
|
||||
}
|
||||
}
|
||||
|
||||
pub fn viscosity_at(&self, temp: f64) -> Option<f64> {
|
||||
if self.viscosities.is_empty() {
|
||||
return None;
|
||||
}
|
||||
let i_viscosity = self
|
||||
.viscosities
|
||||
.iter()
|
||||
.rposition(|viscosity| temp > viscosity.temp_range.0)
|
||||
.unwrap_or(0);
|
||||
Some(self.viscosities[i_viscosity].compute(temp))
|
||||
}
|
||||
|
||||
pub fn conductivity_at(&self, temp: f64) -> Option<f64> {
|
||||
if self.conductivities.is_empty() {
|
||||
return None;
|
||||
}
|
||||
let i_conductivity = self
|
||||
.conductivities
|
||||
.iter()
|
||||
.rposition(|conductivity| temp > conductivity.temp_range.0)
|
||||
.unwrap_or(0);
|
||||
Some(self.conductivities[i_conductivity].compute(temp))
|
||||
}
|
||||
}
|
||||
|
||||
pub struct TransportFit {
|
||||
|
|
@ -11,3 +49,148 @@ pub struct TransportFit {
|
|||
pub c: f64,
|
||||
pub d: f64,
|
||||
}
|
||||
|
||||
impl TransportFit {
|
||||
pub fn compute(&self, temp: f64) -> f64 {
|
||||
let inv_temp = 1.0 / temp;
|
||||
self.a * temp.ln() + self.b * inv_temp + self.c * inv_temp * inv_temp + self.d
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod test {
|
||||
use crate::{
|
||||
assert_delta,
|
||||
properties::transport_fit::{SpeciesTransportData, TransportFit},
|
||||
};
|
||||
|
||||
#[test]
|
||||
fn test_transport_fit_compute() {
|
||||
let fit = TransportFit {
|
||||
temp_range: (1.0, 2.0),
|
||||
a: 10.0,
|
||||
b: 20.0,
|
||||
c: 30.0,
|
||||
d: 40.0,
|
||||
};
|
||||
|
||||
assert_delta!(fit.compute(4.0), 60.73794361119891, 1e-12);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_calc_transport_properties() {
|
||||
let viscosities = vec![
|
||||
TransportFit {
|
||||
temp_range: (1.0, 100.0),
|
||||
a: 10.0,
|
||||
b: 20.0,
|
||||
c: 30.0,
|
||||
d: 40.0,
|
||||
},
|
||||
TransportFit {
|
||||
temp_range: (100.0, 200.0),
|
||||
a: 1.0,
|
||||
b: 2.0,
|
||||
c: 3.0,
|
||||
d: 4.0,
|
||||
},
|
||||
];
|
||||
|
||||
let conductivities = vec![
|
||||
TransportFit {
|
||||
temp_range: (1.0, 100.0),
|
||||
a: 10.0,
|
||||
b: 20.0,
|
||||
c: 30.0,
|
||||
d: 40.0,
|
||||
},
|
||||
TransportFit {
|
||||
temp_range: (100.0, 200.0),
|
||||
a: 1.0,
|
||||
b: 2.0,
|
||||
c: 3.0,
|
||||
d: 4.0,
|
||||
},
|
||||
];
|
||||
|
||||
let data = SpeciesTransportData {
|
||||
viscosities,
|
||||
conductivities,
|
||||
name: "".to_string(),
|
||||
};
|
||||
|
||||
assert_delta!(
|
||||
data.conductivity_at(0.5),
|
||||
data.conductivities[0].compute(0.5),
|
||||
1e-12
|
||||
);
|
||||
assert_delta!(
|
||||
data.conductivity_at(1.0),
|
||||
data.conductivities[0].compute(1.0),
|
||||
1e-12
|
||||
);
|
||||
assert_delta!(
|
||||
data.conductivity_at(50.0),
|
||||
data.conductivities[0].compute(50.0),
|
||||
1e-12
|
||||
);
|
||||
assert_delta!(
|
||||
data.conductivity_at(100.0),
|
||||
data.conductivities[0].compute(100.0),
|
||||
1e-12
|
||||
);
|
||||
assert_delta!(
|
||||
data.conductivity_at(100.0 + 1e-12),
|
||||
data.conductivities[1].compute(100.0 + 1e-12),
|
||||
1e-12
|
||||
);
|
||||
|
||||
assert_delta!(
|
||||
data.conductivity_at(200.0),
|
||||
data.conductivities[1].compute(200.0),
|
||||
1e-12
|
||||
);
|
||||
assert_delta!(
|
||||
data.conductivity_at(500.0),
|
||||
data.conductivities[1].compute(500.0),
|
||||
1e-12
|
||||
);
|
||||
|
||||
assert_delta!(
|
||||
data.viscosity_at(0.5),
|
||||
data.viscosities[0].compute(0.5),
|
||||
1e-12
|
||||
);
|
||||
assert_delta!(
|
||||
data.viscosity_at(1.0),
|
||||
data.viscosities[0].compute(1.0),
|
||||
1e-12
|
||||
);
|
||||
assert_delta!(
|
||||
data.viscosity_at(50.0),
|
||||
data.viscosities[0].compute(50.0),
|
||||
1e-12
|
||||
);
|
||||
assert_delta!(
|
||||
data.viscosity_at(100.0),
|
||||
data.viscosities[0].compute(100.0),
|
||||
1e-12
|
||||
);
|
||||
assert_delta!(
|
||||
data.viscosity_at(100.0 + 1e-12),
|
||||
data.viscosities[1].compute(100.0 + 1e-12),
|
||||
1e-12
|
||||
);
|
||||
|
||||
assert_delta!(
|
||||
data.viscosity_at(200.0),
|
||||
data.viscosities[1].compute(200.0),
|
||||
1e-12
|
||||
);
|
||||
assert_delta!(
|
||||
data.viscosity_at(500.0),
|
||||
data.viscosities[1].compute(500.0),
|
||||
1e-12
|
||||
);
|
||||
}
|
||||
}
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue