fastsim_core/
gas_properties.rs

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

fastsim_core/
gas_properties.rs
