fastsim_core/
gas_properties.rs

1use crate::imports::*;
2
3lazy_static! {
4    /// room temperature
5    pub static ref TE_STD_AIR: si::Temperature = (22. + 273.15) * uc::KELVIN;
6    /// pressure of air at 180 m and 22 C
7    pub static ref STD_PRESSURE_AIR: si::Pressure = 99_346.3 * uc::PASCAL;
8    /// density of air at 180 m ASL and 22 C
9    pub static ref STD_DENSITY_AIR: si::MassDensity = 1.172 * uc::KGPM3;
10    /// ideal gas constant for air
11    pub static ref R_AIR: si::SpecificHeatCapacity = 287.0 * uc::J_PER_KG_K;
12    /// standard elevation above sea level
13    pub static ref H_STD: si::Length = 180.0 * uc::M;
14}
15
16#[serde_api]
17#[derive(Deserialize, Serialize, Debug, Clone, PartialEq)]
18#[non_exhaustive]
19#[cfg_attr(feature = "pyo3", pyclass(module = "fastsim", subclass, eq))]
20#[serde(deny_unknown_fields)]
21pub struct Air {}
22impl Init for Air {}
23impl SerdeAPI for Air {}
24
25#[named_struct_pyo3_api]
26impl Air {
27    #[new]
28    fn __new__() -> Self {
29        Self {}
30    }
31    /// Returns density of air \[kg/m^3\]
32    /// Source: <https://www.grc.nasa.gov/WWW/K-12/rocket/atmosmet.html>  
33    ///
34    /// # Equations used
35    /// T = 15.04 - .00649 * h  
36    /// p = 101.29 * [(T + 273.1)/288.08]^5.256  
37    ///
38    /// # Arguments  
39    /// * `te_air_deg_c` - optional ambient temperature [°C] of air, defaults to 22 C
40    /// * `h_m` - optional elevation \[m\] above sea level, defaults to 180 m
41    #[staticmethod]
42    #[pyo3(name = "get_density")]
43    #[pyo3(signature = (te_air_deg_c=None, h_m=None))]
44    pub fn get_density_py(te_air_deg_c: Option<f64>, h_m: Option<f64>) -> f64 {
45        Self::get_density(
46            te_air_deg_c.map(|te_air_deg_c| (te_air_deg_c + 273.15) * uc::KELVIN),
47            h_m.map(|h_m| h_m * uc::M),
48        )
49        .get::<si::kilogram_per_cubic_meter>()
50    }
51
52    /// Returns thermal conductivity [W/(m*K)] of air
53    /// # Arguments
54    /// - `te_air`: temperature [°C] of air
55    #[pyo3(name = "get_therm_cond")]
56    #[staticmethod]
57    pub fn get_therm_cond_py(te_air: f64) -> anyhow::Result<f64> {
58        Ok(
59            Self::get_therm_cond((te_air + uc::CELSIUS_TO_KELVIN) * uc::KELVIN)?
60                .get::<si::watt_per_meter_kelvin>(),
61        )
62    }
63
64    /// Returns constant pressure specific heat [J/(kg*K)] of air
65    /// # Arguments
66    /// - `te_air`: temperature [°C] of air
67    #[pyo3(name = "get_specific_heat_cp")]
68    #[staticmethod]
69    pub fn get_specific_heat_cp_py(te_air: f64) -> anyhow::Result<f64> {
70        Ok(
71            Self::get_specific_heat_cp((te_air + uc::CELSIUS_TO_KELVIN) * uc::KELVIN)?
72                .get::<si::joule_per_kilogram_kelvin>(),
73        )
74    }
75
76    /// Returns specific enthalpy [J/kg] of air  
77    /// # Arguments  
78    /// - `te_air`: temperature [°C] of air
79    #[pyo3(name = "get_specific_enthalpy")]
80    #[staticmethod]
81    pub fn get_specific_enthalpy_py(te_air: f64) -> anyhow::Result<f64> {
82        Ok(
83            Self::get_specific_enthalpy((te_air - uc::CELSIUS_TO_KELVIN) * uc::KELVIN)?
84                .get::<si::joule_per_kilogram>(),
85        )
86    }
87
88    /// Returns specific energy [J/kg] of air  
89    /// # Arguments  
90    /// - `te_air`: temperature [°C] of air
91    #[pyo3(name = "get_specific_energy")]
92    #[staticmethod]
93    pub fn get_specific_energy_py(te_air: f64) -> anyhow::Result<f64> {
94        Ok(
95            Self::get_specific_energy((te_air - uc::CELSIUS_TO_KELVIN) * uc::KELVIN)?
96                .get::<si::joule_per_kilogram>(),
97        )
98    }
99
100    /// Returns thermal Prandtl number of air
101    /// # Arguments
102    /// - `te_air`: temperature [°C] of air     
103    #[pyo3(name = "get_pr")]
104    #[staticmethod]
105    pub fn get_pr_py(te_air: f64) -> anyhow::Result<f64> {
106        Ok(Self::get_pr((te_air - uc::CELSIUS_TO_KELVIN) * uc::KELVIN)?.get::<si::ratio>())
107    }
108
109    /// Returns dynamic viscosity \[Pa*s\] of air
110    /// # Arguments
111    /// te_air: temperature [°C] of air
112    #[pyo3(name = "get_dyn_visc")]
113    #[staticmethod]
114    pub fn get_dyn_visc_py(te_air: f64) -> anyhow::Result<f64> {
115        Ok(
116            Self::get_dyn_visc((te_air - uc::CELSIUS_TO_KELVIN) * uc::KELVIN)?
117                .get::<si::pascal_second>(),
118        )
119    }
120
121    /// Returns temperature [°C] of air
122    /// # Arguments
123    /// - `h`: specific enthalpy of air \[J/kg\]
124    #[pyo3(name = "get_te_from_h")]
125    #[staticmethod]
126    pub fn get_te_from_h_py(h: f64) -> anyhow::Result<f64> {
127        Ok(Self::get_te_from_h(h * uc::J_PER_KG)?.get::<si::degree_celsius>())
128    }
129
130    /// Returns temperature [°C] of air
131    /// # Arguments
132    /// - `u`: specific energy of air \[J/kg\]
133    #[pyo3(name = "get_te_from_u")]
134    #[staticmethod]
135    pub fn get_te_from_u_py(u: f64) -> anyhow::Result<f64> {
136        Ok(Self::get_te_from_u(u * uc::J_PER_KG)?.get::<si::degree_celsius>())
137    }
138}
139
140impl Air {
141    /// Returns density of air
142    /// Source: <https://www.grc.nasa.gov/WWW/K-12/rocket/atmosmet.html>  
143    /// Note that if `None` is passed for either argument, function evaluation should be faster
144    ///
145    /// # Equations used
146    /// - T = 15.04 - 0.00649 * h  
147    /// - p = 101.29 * ((T + 273.1) / 288.08) ^ 5.256  
148    ///
149    /// # Arguments  
150    /// * `te_air` - ambient temperature of air, defaults to 22 C
151    /// * `h` - elevation above sea level, defaults to 180 m
152    pub fn get_density(te_air: Option<si::Temperature>, h: Option<si::Length>) -> si::MassDensity {
153        let std_pressure_at_elev = |h: si::Length| -> si::Pressure {
154            let std_temp_at_elev = (15.04 - 0.00649 * h.get::<si::meter>() + 273.15) * uc::KELVIN;
155            (101.29e3 * uc::PASCAL)
156                * ((std_temp_at_elev / (288.08 * uc::KELVIN))
157                    .get::<si::ratio>()
158                    .powf(5.256))
159        };
160        match (h, te_air) {
161            (None, None) => *STD_DENSITY_AIR,
162            (None, Some(te_air)) => *STD_PRESSURE_AIR / *R_AIR / te_air,
163            (Some(h_val), None) => std_pressure_at_elev(h_val) / *R_AIR / *TE_STD_AIR,
164            (Some(h_val), Some(te_air)) => std_pressure_at_elev(h_val) / *R_AIR / te_air,
165        }
166    }
167
168    /// Returns thermal conductivity of air
169    /// # Arguments
170    /// - `te_air`: temperature of air
171    pub fn get_therm_cond(te_air: si::Temperature) -> anyhow::Result<si::ThermalConductivity> {
172        Ok(
173            asp::THERMAL_CONDUCTIVITY_INTERP.interpolate(&[te_air.get::<si::degree_celsius>()])?
174                * uc::WATT_PER_METER_KELVIN,
175        )
176    }
177
178    /// Returns constant pressure specific heat of air
179    /// # Arguments
180    /// - `te_air`: temperature of air
181    pub fn get_specific_heat_cp(
182        te_air: si::Temperature,
183    ) -> anyhow::Result<si::SpecificHeatCapacity> {
184        Ok(asp::C_P_INTERP.interpolate(&[te_air.get::<si::degree_celsius>()])? * uc::J_PER_KG_K)
185    }
186
187    /// Returns specific enthalpy of air  
188    /// # Arguments  
189    /// - `te_air`: temperature of air
190    pub fn get_specific_enthalpy(te_air: si::Temperature) -> anyhow::Result<si::SpecificEnergy> {
191        Ok(asp::ENTHALPY_INTERP.interpolate(&[te_air.get::<si::degree_celsius>()])? * uc::J_PER_KG)
192    }
193
194    /// Returns specific energy of air  
195    /// # Arguments  
196    /// - `te_air`: temperature of air
197    pub fn get_specific_energy(te_air: si::Temperature) -> anyhow::Result<si::SpecificEnergy> {
198        Ok(asp::ENERGY_INTERP.interpolate(&[te_air.get::<si::degree_celsius>()])? * uc::J_PER_KG)
199    }
200
201    /// Returns thermal Prandtl number of air
202    /// # Arguments
203    /// - `te_air`: temperature of air     
204    pub fn get_pr(te_air: si::Temperature) -> anyhow::Result<si::Ratio> {
205        Ok(asp::PRANDTL_INTERP.interpolate(&[te_air.get::<si::degree_celsius>()])? * uc::R)
206    }
207
208    /// Returns dynamic viscosity \[Pa*s\] of air
209    /// # Arguments
210    /// te_air: temperature of air
211    pub fn get_dyn_visc(te_air: si::Temperature) -> anyhow::Result<si::DynamicViscosity> {
212        Ok(
213            asp::DYN_VISC_INTERP.interpolate(&[te_air.get::<si::degree_celsius>()])?
214                * uc::PASCAL_SECOND,
215        )
216    }
217
218    /// Returns temperature of air
219    /// # Arguments
220    /// `h`: specific enthalpy of air \[J/kg\]
221    pub fn get_te_from_h(h: si::SpecificEnergy) -> anyhow::Result<si::Temperature> {
222        Ok(asp::TEMP_FROM_ENTHALPY.interpolate(&[h.get::<si::joule_per_kilogram>()])? * uc::KELVIN)
223    }
224
225    /// Returns temperature of air
226    /// # Arguments
227    /// `u`: specific energy of air \[J/kg\]
228    pub fn get_te_from_u(u: si::SpecificEnergy) -> anyhow::Result<si::Temperature> {
229        Ok(asp::TEMP_FROM_ENERGY.interpolate(&[u.get::<si::joule_per_kilogram>()])? * uc::KELVIN)
230    }
231}
232
233use air_static_props as asp;
234
235/// Air fluid properties for calculations.  
236///
237/// Values obtained via (in Python, after running `pip install CoolProp`):
238/// ```python
239/// from CoolProp.CoolProp import PropsSI
240/// import numpy as np
241/// import pandas as pd
242/// T_degC = np.logspace(1, np.log10(5e3 + 70), 25) - 70
243/// T = T_degC + 273.15
244/// prop_dict = {
245///     'T [°C]': T_degC,
246///     'h [J/kg]': [0] * len(T),
247///     'u [J/kg]': [0] * len(T),
248///     'k [W/(m*K)]': [0] * len(T),
249///     'rho [kg/m^3]': [0] * len(T),
250///     'c_p [J/(kg*K)]': [0] * len(T),
251///     'mu [Pa*s]': [0] * len(T),
252/// }
253///
254/// species = "Air"
255///
256/// for i, _ in enumerate(T_degC):
257///     prop_dict['h [J/kg]'][i] = f"{PropsSI('H', 'P', 101325, 'T', T[i], species):.5g}" # specific enthalpy [J/(kg*K)]
258///     prop_dict['u [J/kg]'][i] = f"{PropsSI('U', 'P', 101325, 'T', T[i], species):.5g}" # specific enthalpy [J/(kg*K)]
259///     prop_dict['k [W/(m*K)]'][i] = f"{PropsSI('L', 'P', 101325, 'T', T[i], species):.5g}" # thermal conductivity [W/(m*K)]
260///     prop_dict['rho [kg/m^3]'][i] = f"{PropsSI('D', 'P', 101325, 'T', T[i], species):.5g}" # density [kg/m^3]
261///     prop_dict['c_p [J/(kg*K)]'][i] = f"{PropsSI('C', 'P', 101325, 'T', T[i], species):.5g}" # density [kg/m^3]
262///     prop_dict['mu [Pa*s]'][i] = f"{PropsSI('V', 'P', 101325, 'T', T[i], species):.5g}" # viscosity [Pa*s]
263///
264/// prop_df = pd.DataFrame(data=prop_dict)
265/// pd.set_option('display.float_format', lambda x: '%.3g' % x)
266/// prop_df = prop_df.apply(np.float64)
267/// ```
268mod air_static_props {
269    use super::*;
270    lazy_static! {
271        /// Array of temperatures at which properties are evaluated
272        static ref TEMPERATURE_DEG_C_VALUES: Vec<f64> = vec![
273            -60.,
274             -57.03690616,
275            -53.1958198,
276            -48.21658352,
277            -41.7619528,
278            -33.39475442,
279            -22.54827664,
280            -8.48788571,
281            9.73873099,
282            33.36606527,
283            63.99440042,
284            103.69819869,
285            155.16660498,
286            221.88558305,
287            308.37402042,
288            420.48979341,
289            565.82652205,
290            754.22788725,
291            998.45434496,
292            1315.04739396,
293            1725.44993435,
294            2257.45859876,
295            2947.10642291,
296            3841.10336915,
297            5000.
298        ];
299        pub static ref TEMP_FROM_ENTHALPY: Interpolator = Interpolator::new_1d(
300            ENTHALPY_VALUES.iter().map(|x| x.get::<si::joule_per_kilogram>()).collect::<Vec<f64>>(),
301            TEMPERATURE_DEG_C_VALUES.clone(),
302            Strategy::Linear,
303            Extrapolate::Error
304        ).unwrap_or_else(|_| panic!("Failed to construct gas properties vec"));
305        pub static ref TEMP_FROM_ENERGY: Interpolator = Interpolator::new_1d(
306            ENERGY_VALUES.iter().map(|x| x.get::<si::joule_per_kilogram>()).collect::<Vec<f64>>(),
307            TEMPERATURE_DEG_C_VALUES.clone(),
308            Strategy::Linear,
309            Extrapolate::Error
310        ).unwrap_or_else(|_| panic!("Failed to construct gas properties vec"));
311        /// Thermal conductivity values of air corresponding to temperature values
312        static ref THERMAL_CONDUCTIVITY_VALUES: Vec<si::ThermalConductivity> = [
313            0.019597,
314            0.019841,
315            0.020156,
316            0.020561,
317            0.021083,
318            0.021753,
319            0.022612,
320            0.023708,
321            0.025102,
322            0.026867,
323            0.02909,
324            0.031875,
325            0.035342,
326            0.039633,
327            0.044917,
328            0.051398,
329            0.059334,
330            0.069059,
331            0.081025,
332            0.095855,
333            0.11442,
334            0.13797,
335            0.16828,
336            0.20795,
337            0.26081
338        ]
339        .iter()
340        .map(|x| *x * uc::WATT_PER_METER_KELVIN)
341        .collect();
342        pub static ref THERMAL_CONDUCTIVITY_INTERP: Interpolator = Interpolator::new_1d(
343            TEMPERATURE_DEG_C_VALUES.clone(),
344            THERMAL_CONDUCTIVITY_VALUES.iter().map(|x| x.get::<si::watt_per_meter_degree_celsius>()).collect::<Vec<f64>>(),
345            Strategy::Linear,
346            Extrapolate::Error
347        ).unwrap_or_else(|_| panic!("Failed to construct gas properties vec"));
348        /// Specific heat values of air corresponding to temperature values
349        static ref C_P_VALUES: Vec<si::SpecificHeatCapacity> = [
350            1006.2,
351            1006.1,
352            1006.,
353            1005.9,
354            1005.7,
355            1005.6,
356            1005.5,
357            1005.6,
358            1005.9,
359            1006.6,
360            1008.3,
361            1011.6,
362            1017.9,
363            1028.9,
364            1047.,
365            1073.4,
366            1107.6,
367            1146.1,
368            1184.5,
369            1219.5,
370            1250.1,
371            1277.1,
372            1301.7,
373            1324.5,
374            1347.
375        ]
376        .iter()
377        .map(|x| *x * uc::J_PER_KG_K)
378        .collect();
379        pub static ref C_P_INTERP: Interpolator = Interpolator::new_1d(
380            TEMPERATURE_DEG_C_VALUES.clone(),
381            C_P_VALUES.iter().map(|x| x.get::<si::joule_per_kilogram_kelvin>()).collect::<Vec<f64>>(),
382            Strategy::Linear,
383            Extrapolate::Error
384        ).unwrap_or_else(|_| panic!("Failed to construct gas properties vec"));
385        static ref ENTHALPY_VALUES: Vec<si::SpecificEnergy> = [
386            338940.,
387            341930.,
388            345790.,
389            350800.,
390            357290.,
391            365710.,
392            376610.,
393            390750.,
394            409080.,
395            432860.,
396            463710.,
397            503800.,
398            556020.,
399            624280.,
400            714030.,
401            832880.,
402            991400.,
403            1203800.,
404            1488700.,
405            1869600.,
406            2376700.,
407            3049400.,
408            3939100.,
409            5113600.,
410            6662000.
411        ]
412        .iter()
413        .map(|x| *x * uc::J_PER_KG)
414        .collect();
415        pub static ref ENTHALPY_INTERP: Interpolator = Interpolator::new_1d(
416            TEMPERATURE_DEG_C_VALUES.clone(),
417            ENTHALPY_VALUES.iter().map(|x| x.get::<si::joule_per_kilogram>()).collect::<Vec<f64>>(),
418            Strategy::Linear,
419            Extrapolate::Error
420        ).unwrap_or_else(|_| panic!("Failed to construct gas properties vec"));
421        pub static ref ENERGY_VALUES: Vec<si::SpecificEnergy> = [
422            277880.,
423            280000.,
424            282760.,
425            286330.,
426            290960.,
427            296960.,
428            304750.,
429            314840.,
430            327920.,
431            344890.,
432            366940.,
433            395620.,
434            433040.,
435            482140.,
436            547050.,
437            633700.,
438            750490.,
439            908830.,
440            1123600.,
441            1413600.,
442            1802900.,
443            2322900.,
444            3014700.,
445            3932500.,
446            5148300.
447        ]
448        .iter()
449        .map(|x| *x * uc::J_PER_KG)
450        .collect();
451        pub static ref ENERGY_INTERP: Interpolator = Interpolator::new_1d(
452            TEMPERATURE_DEG_C_VALUES.clone(),
453            ENERGY_VALUES.iter().map(|x| x.get::<si::joule_per_kilogram>()).collect::<Vec<f64>>(),
454            Strategy::Linear,
455            Extrapolate::Error
456        ).unwrap_or_else(|_| panic!("Failed to construct gas properties vec"));
457        static ref DYN_VISCOSITY_VALUES: Vec<si::DynamicViscosity> = [
458            1.4067e-05,
459            1.4230e-05,
460            1.4440e-05,
461            1.4711e-05,
462            1.5058e-05,
463            1.5502e-05,
464            1.6069e-05,
465            1.6791e-05,
466            1.7703e-05,
467            1.8850e-05,
468            2.0283e-05,
469            2.2058e-05,
470            2.4240e-05,
471            2.6899e-05,
472            3.0112e-05,
473            3.3966e-05,
474            3.8567e-05,
475            4.4049e-05,
476            5.0595e-05,
477            5.8464e-05,
478            6.8036e-05,
479            7.9878e-05,
480            9.4840e-05,
481            1.1423e-04,
482            1.4006e-04
483        ]
484        .iter()
485        .map(|x| *x * uc::PASCAL_SECOND)
486        .collect();
487        pub static ref DYN_VISC_INTERP: Interpolator = Interpolator::new_1d(
488            TEMPERATURE_DEG_C_VALUES.clone(),
489            DYN_VISCOSITY_VALUES.iter().map(|x| x.get::<si::pascal_second>()).collect::<Vec<f64>>(),
490            Strategy::Linear,
491            Extrapolate::Error
492        ).unwrap_or_else(|_| panic!("Failed to construct gas properties vec"));
493        static ref PRANDTL_VALUES: Vec<si::Ratio> = DYN_VISCOSITY_VALUES
494            .iter()
495            .zip(C_P_VALUES.iter())
496            .zip(THERMAL_CONDUCTIVITY_VALUES.iter())
497            .map(|((mu, c_p), k)| -> si::Ratio {*mu * *c_p / *k})
498            .collect::<Vec<si::Ratio>>();
499        pub static ref PRANDTL_INTERP: Interpolator = Interpolator::new_1d(
500            TEMPERATURE_DEG_C_VALUES.clone(),
501            PRANDTL_VALUES.iter().map(|x| x.get::<si::ratio>()).collect::<Vec<f64>>(),
502            Strategy::Linear,
503            Extrapolate::Error
504        ).unwrap_or_else(|_| panic!("Failed to construct gas properties vec"));
505    }
506}
507
508use octane_static_props as osp;
509
510/// Octane (as a surrogate for gasoline) fluid properties for calculations.  
511///
512/// Values obtained via (in Python, after running `pip install CoolProp`):
513/// ```python
514/// from CoolProp.CoolProp import PropsSI
515/// import numpy as np
516/// import pandas as pd
517/// T_degC = np.logspace(1, np.log10(5e3 + 70), 25) - 50
518/// T = T_degC + 273.15
519/// prop_dict = {
520///     'T [°C]': T_degC,
521///     'h [J/kg]': [0] * len(T),
522///     'u [J/kg]': [0] * len(T),
523///     'k [W/(m*K)]': [0] * len(T),
524///     'rho [kg/m^3]': [0] * len(T),
525///     'c_p [J/(kg*K)]': [0] * len(T),
526///     'mu [Pa*s]': [0] * len(T),
527/// }
528///
529/// species = "Octane"
530///
531/// for i, _ in enumerate(T_degC):
532///     prop_dict['h [J/kg]'][i] = f"{PropsSI('H', 'P', 101325, 'T', T[i], species):.5g}" # specific enthalpy [J/(kg*K)]
533///     prop_dict['u [J/kg]'][i] = f"{PropsSI('U', 'P', 101325, 'T', T[i], species):.5g}" # specific enthalpy [J/(kg*K)]
534///     prop_dict['k [W/(m*K)]'][i] = f"{PropsSI('L', 'P', 101325, 'T', T[i], species):.5g}" # thermal conductivity [W/(m*K)]
535///     prop_dict['rho [kg/m^3]'][i] = f"{PropsSI('D', 'P', 101325, 'T', T[i], species):.5g}" # density [kg/m^3]
536///     prop_dict['c_p [J/(kg*K)]'][i] = f"{PropsSI('C', 'P', 101325, 'T', T[i], species):.5g}" # density [kg/m^3]
537///     prop_dict['mu [Pa*s]'][i] = f"{PropsSI('V', 'P', 101325, 'T', T[i], species):.5g}" # viscosity [Pa*s]
538///
539/// prop_df = pd.DataFrame(data=prop_dict)
540/// pd.set_option('display.float_format', lambda x: '%.3g' % x)
541/// prop_df = prop_df.apply(np.float64)
542/// ```
543mod octane_static_props {
544    use super::*;
545    lazy_static! {
546        /// Array of temperatures at which properties are evaluated
547        static ref TEMPERATURE_DEG_C_VALUES: Vec<f64> = vec![
548            -4.00000000e+01,
549             -3.70369062e+01,
550            -3.31958198e+01,
551            -2.82165835e+01,
552            -2.17619528e+01,
553            -1.33947544e+01,
554            -2.54827664e+00,
555            1.15121143e+01,
556            2.97387310e+01,
557            5.33660653e+01,
558            8.39944004e+01,
559            1.23698199e+02,
560            1.75166605e+02,
561            2.41885583e+02,
562            3.28374020e+02,
563            4.40489793e+02,
564            5.85826522e+02,
565            7.74227887e+02,
566            1.01845434e+03,
567            1.33504739e+03,
568            1.74544993e+03,
569            2.27745860e+03,
570            2.96710642e+03,
571            3.86110337e+03,
572            5.02000000e+03
573        ];
574        pub static ref TEMP_FROM_ENERGY: Interpolator = Interpolator::new_1d(
575           ENERGY_VALUES.iter().map(|x| x.get::<si::joule_per_kilogram>()).collect::<Vec<f64>>(),
576           TEMPERATURE_DEG_C_VALUES.clone(),
577           Strategy::Linear,
578           Extrapolate::Error
579        ).unwrap_or_else(|_| panic!("Failed to construct gas properties vec"));
580        pub static ref ENERGY_VALUES: Vec<si::SpecificEnergy> = [
581            -3.8247e+05,
582            -3.7645e+05,
583            -3.6862e+05,
584            -3.5841e+05,
585            -3.4507e+05,
586            -3.2760e+05,
587            -3.0464e+05,
588            -2.7432e+05,
589            -2.3400e+05,
590            -1.7991e+05,
591            -1.0649e+05,
592            -5.3074e+03,
593            3.8083e+05,
594            5.3958e+05,
595            7.6926e+05,
596            1.1024e+06,
597            1.5836e+06,
598            2.2729e+06,
599            3.2470e+06,
600            4.6015e+06,
601            6.4541e+06,
602            8.9500e+06,
603            1.2272e+07,
604            1.6654e+07,
605            2.2399e+07
606        ]
607        .iter()
608        .map(|x| *x * uc::J_PER_KG)
609        .collect();
610        pub static ref ENERGY_INTERP: Interpolator = Interpolator::new_1d(
611           TEMPERATURE_DEG_C_VALUES.clone(),
612           ENERGY_VALUES.iter().map(|x| x.get::<si::joule_per_kilogram>()).collect::<Vec<f64>>(),
613           Strategy::Linear,
614           Extrapolate::Error
615        ).unwrap_or_else(|_| panic!("Failed to construct gas properties vec"));
616    }
617}
618
619#[serde_api]
620#[derive(Deserialize, Serialize, Debug, Clone, PartialEq, StateMethods)]
621#[cfg_attr(feature = "pyo3", pyclass(module = "fastsim", subclass, eq))]
622#[serde(deny_unknown_fields)]
623pub struct Octane {}
624impl Init for Octane {}
625impl SerdeAPI for Octane {}
626
627#[named_struct_pyo3_api]
628impl Octane {
629    /// Returns specific energy [J/kg] of octane  
630    /// # Arguments  
631    /// - `te_octane`: temperature [°C] of octane
632    #[pyo3(name = "get_specific_energy")]
633    #[staticmethod]
634    pub fn get_specific_energy_py(te_octane: f64) -> anyhow::Result<f64> {
635        Ok(
636            Self::get_specific_energy((te_octane - uc::CELSIUS_TO_KELVIN) * uc::KELVIN)?
637                .get::<si::joule_per_kilogram>(),
638        )
639    }
640
641    /// Returns temperature [°C] of octane
642    /// # Arguments
643    /// - `u`: specific energy of octane \[J/kg\]
644    #[pyo3(name = "get_te_from_u")]
645    #[staticmethod]
646    pub fn get_te_from_u_py(u: f64) -> anyhow::Result<f64> {
647        Ok(Self::get_te_from_u(u * uc::J_PER_KG)?.get::<si::degree_celsius>())
648    }
649}
650
651impl Octane {
652    /// Returns specific energy of octane  
653    /// # Arguments  
654    /// - `te_octane`: temperature of octane
655    pub fn get_specific_energy(te_octane: si::Temperature) -> anyhow::Result<si::SpecificEnergy> {
656        Ok(
657            osp::ENERGY_INTERP.interpolate(&[te_octane.get::<si::degree_celsius>()])?
658                * uc::J_PER_KG,
659        )
660    }
661
662    /// Returns temperature of octane
663    /// # Arguments
664    /// `u`: specific energy of octane \[J/kg\]
665    pub fn get_te_from_u(u: si::SpecificEnergy) -> anyhow::Result<si::Temperature> {
666        Ok(osp::TEMP_FROM_ENERGY.interpolate(&[u.get::<si::joule_per_kilogram>()])? * uc::KELVIN)
667    }
668}
669
670/// Given Reynolds number `re`, return C and m to calculate Nusselt number for
671/// sphere, from Incropera's Intro to Heat Transfer, 5th Ed., eq. 7.44
672pub fn get_sphere_conv_params(re: f64) -> (f64, f64) {
673    let (c, m) = if re < 4.0 {
674        (0.989, 0.330)
675    } else if re < 40.0 {
676        (0.911, 0.385)
677    } else if re < 4e3 {
678        (0.683, 0.466)
679    } else if re < 40e3 {
680        (0.193, 0.618)
681    } else {
682        (0.027, 0.805)
683    };
684    (c, m)
685}
fastsim_core/gas_properties.rs

fastsim_core/
gas_properties.rs
