fastsim_core/
thermal.rs

1//! Module for simulating thermal behavior of powertrains
2
3use crate::proc_macros::{add_pyo3_api, HistoryVec};
4
5use crate::air::AirProperties;
6use crate::cycle;
7use crate::imports::*;
8#[cfg(feature = "pyo3")]
9use crate::pyo3imports::*;
10use crate::simdrive;
11use crate::vehicle;
12use crate::vehicle_thermal::*;
13
14#[add_pyo3_api(
15    /// method for instantiating SimDriveHot
16    #[new]
17    #[pyo3(signature = (cyc, veh, vehthrm, init_state=None, amb_te_deg_c=None))]
18    pub fn __new__(
19        cyc: cycle::RustCycle,
20        veh: vehicle::RustVehicle,
21        vehthrm: VehicleThermal,
22        init_state: Option<ThermalState>,
23        amb_te_deg_c: Option<Vec<f64>>,
24     ) -> Self {
25        Self::new(cyc, veh, vehthrm, init_state, amb_te_deg_c.map(Array1::from))
26    }
27
28    #[pyo3(name = "gap_to_lead_vehicle_m")]
29    /// Provides the gap-with lead vehicle from start to finish
30    pub fn gap_to_lead_vehicle_m_py(&self) -> anyhow::Result<Vec<f64>> {
31        Ok(self.gap_to_lead_vehicle_m().to_vec())
32    }
33     #[pyo3(name = "sim_drive")]
34     #[pyo3(signature = (init_soc=None, aux_in_kw_override=None))]
35    /// Initialize and run sim_drive_walk as appropriate for vehicle attribute vehPtType.
36    /// Arguments
37    /// ------------
38    /// init_soc: initial SOC for electrified vehicles.
39    /// aux_in_kw: aux_in_kw override.  Array of same length as cyc.time_s.
40    ///     Default of None causes veh.aux_kw to be used.
41    pub fn sim_drive_py(
42        &mut self,
43        init_soc: Option<f64>,
44        aux_in_kw_override: Option<Vec<f64>>,
45    ) -> anyhow::Result<()> {
46        let aux_in_kw_override = aux_in_kw_override.map(Array1::from);
47        self.sim_drive(init_soc, aux_in_kw_override)
48    }
49
50    #[pyo3(signature = (init_soc, aux_in_kw_override=None))]
51    /// Receives second-by-second cycle information, vehicle properties,
52    /// and an initial state of charge and runs sim_drive_step to perform a
53    /// backward facing powertrain simulation. Method 'sim_drive' runs this
54    /// iteratively to achieve correct SOC initial and final conditions, as
55    /// needed.
56    ///
57    /// Arguments
58    /// ------------
59    /// init_soc (optional): initial battery state-of-charge (SOC) for electrified vehicles
60    /// aux_in_kw: aux_in_kw override.  Array of same length as cyc.time_s.
61    ///         None causes veh.aux_kw to be used.
62    pub fn sim_drive_walk(
63        &mut self,
64        init_soc: f64,
65        aux_in_kw_override: Option<Vec<f64>>,
66    ) {
67        let aux_in_kw_override = aux_in_kw_override.map(Array1::from);
68        self.walk(init_soc, aux_in_kw_override);
69    }
70
71    #[pyo3(name = "init_for_step")]
72    #[pyo3(signature = (init_soc, aux_in_kw_override=None))]
73    /// This is a specialty method which should be called prior to using
74    /// sim_drive_step in a loop.
75    /// Arguments
76    /// ------------
77    /// init_soc: initial battery state-of-charge (SOC) for electrified vehicles
78    /// aux_in_kw: aux_in_kw override.  Array of same length as cyc.time_s.
79    ///         Default of None causes veh.aux_kw to be used.
80    pub fn init_for_step_py(
81        &mut self,
82        init_soc:f64,
83        aux_in_kw_override: Option<Vec<f64>>
84    ) {
85        let aux_in_kw_override = aux_in_kw_override.map(Array1::from);
86        self.init_for_step(init_soc, aux_in_kw_override);
87    }
88
89    /// Step through 1 time step.
90    pub fn sim_drive_step(&mut self) -> anyhow::Result<()> {
91        self.step()
92    }
93    #[pyo3(name = "solve_step")]
94    /// Perform all the calculations to solve 1 time step.
95    pub fn solve_step_py(&mut self, i: usize) -> anyhow::Result<()> {
96        self.solve_step(i)
97    }
98
99    #[pyo3(name = "set_misc_calcs")]
100    /// Sets misc. calculations at time step 'i'
101    /// Arguments:
102    /// ----------
103    /// i: index of time step
104    pub fn set_misc_calcs_py(&mut self, i: usize) {
105        self.set_misc_calcs(i);
106    }
107
108    #[pyo3(name = "set_comp_lims")]
109    // Calculate actual speed achieved if vehicle hardware cannot achieve trace speed.
110    // Arguments
111    // ------------
112    // i: index of time step
113    pub fn set_comp_lims_py(&mut self, i: usize) -> anyhow::Result<()> {
114        self.set_comp_lims(i)
115    }
116
117    #[pyo3(name = "set_power_calcs")]
118    /// Calculate power requirements to meet cycle and determine if
119    /// cycle can be met.
120    /// Arguments
121    /// ------------
122    /// i: index of time step
123    pub fn set_power_calcs_py(&mut self, i: usize) -> anyhow::Result<()> {
124        self.set_power_calcs(i)
125    }
126
127    #[pyo3(name = "set_ach_speed")]
128    // Calculate actual speed achieved if vehicle hardware cannot achieve trace speed.
129    // Arguments
130    // ------------
131    // i: index of time step
132    pub fn set_ach_speed_py(&mut self, i: usize) -> anyhow::Result<()> {
133        self.set_ach_speed(i)
134    }
135
136    #[pyo3(name = "set_hybrid_cont_calcs")]
137    /// Hybrid control calculations.
138    /// Arguments
139    /// ------------
140    /// i: index of time step
141    pub fn set_hybrid_cont_calcs_py(&mut self, i: usize) -> anyhow::Result<()> {
142        self.set_hybrid_cont_calcs(i)
143    }
144
145    #[pyo3(name = "set_fc_forced_state")]
146    /// Calculate control variables related to engine on/off state
147    /// Arguments
148    /// ------------
149    /// i: index of time step
150    /// `_py` extension is needed to avoid name collision with getter/setter methods
151    pub fn set_fc_forced_state_py(&mut self, i: usize) -> anyhow::Result<()> {
152        self.set_fc_forced_state_rust(i)
153    }
154
155    #[pyo3(name = "set_hybrid_cont_decisions")]
156    /// Hybrid control decisions.
157    /// Arguments
158    /// ------------
159    /// i: index of time step
160    pub fn set_hybrid_cont_decisions_py(&mut self, i: usize) -> anyhow::Result<()> {
161        self.set_hybrid_cont_decisions(i)
162    }
163
164    #[pyo3(name = "set_fc_power")]
165    /// Sets power consumption values for the current time step.
166    /// Arguments
167    /// ------------
168    /// i: index of time step
169    pub fn set_fc_power_py(&mut self, i: usize) -> anyhow::Result<()> {
170        self.set_fc_power(i)
171    }
172
173    #[pyo3(name = "set_time_dilation")]
174    /// Sets the time dilation for the current step.
175    /// Arguments
176    /// ------------
177    /// i: index of time step
178    pub fn set_time_dilation_py(&mut self, i: usize) -> anyhow::Result<()> {
179        self.set_time_dilation(i)
180    }
181
182    #[pyo3(name = "set_post_scalars")]
183    /// Sets scalar variables that can be calculated after a cycle is run.
184    /// This includes mpgge, various energy metrics, and others
185    pub fn set_post_scalars_py(&mut self) -> anyhow::Result<()> {
186        self.set_post_scalars()
187    }
188)]
189#[derive(Serialize, Deserialize, Clone, PartialEq, Debug)]
190pub struct SimDriveHot {
191    #[api(has_orphaned)]
192    pub sd: simdrive::RustSimDrive,
193    #[api(has_orphaned)]
194    pub vehthrm: VehicleThermal,
195    #[api(skip_get, skip_set)]
196    #[serde(skip)]
197    air: AirProperties,
198    #[api(has_orphaned)]
199    pub state: ThermalState,
200    pub history: ThermalStateHistoryVec,
201    pub hvac_model_history: HVACModelHistoryVec,
202    #[api(skip_get, skip_set)]
203    amb_te_deg_c: Option<Array1<f64>>,
204}
205
206impl SerdeAPI for SimDriveHot {
207    fn init(&mut self) -> anyhow::Result<()> {
208        self.sd.veh.init()?;
209        Ok(())
210    }
211}
212
213impl SimDriveHot {
214    pub fn new(
215        cyc: cycle::RustCycle,
216        veh: vehicle::RustVehicle,
217        vehthrm: VehicleThermal,
218        init_state: Option<ThermalState>,
219        amb_te_deg_c: Option<Array1<f64>>,
220    ) -> Self {
221        let sd = simdrive::RustSimDrive::new(cyc, veh);
222        let air = AirProperties::default();
223        let history = ThermalStateHistoryVec::default();
224
225        let (amb_te_deg_c_arr, state) = match amb_te_deg_c {
226            Some(amb_te_deg_c_arr) => match init_state {
227                Some(state) => {
228                    assert_eq!(state.amb_te_deg_c, amb_te_deg_c_arr[0]);
229                    (Some(amb_te_deg_c_arr), state)
230                }
231                None => {
232                    let state = ThermalState {
233                        amb_te_deg_c: amb_te_deg_c_arr[0],
234                        ..ThermalState::default()
235                    };
236                    (Some(amb_te_deg_c_arr), state)
237                }
238            },
239            None => (
240                None, // 1st return element
241                init_state.unwrap_or_default(),
242            ),
243        };
244
245        Self {
246            sd,
247            vehthrm,
248            air,
249            state,
250            history,
251            hvac_model_history: HVACModelHistoryVec::default(),
252            amb_te_deg_c: amb_te_deg_c_arr,
253        }
254    }
255
256    pub fn gap_to_lead_vehicle_m(&self) -> Array1<f64> {
257        self.sd.gap_to_lead_vehicle_m()
258    }
259
260    pub fn sim_drive(
261        &mut self,
262        init_soc: Option<f64>,
263        aux_in_kw_override: Option<Array1<f64>>,
264    ) -> anyhow::Result<()> {
265        self.sd.hev_sim_count = 0;
266
267        let init_soc = match init_soc {
268            Some(x) => x,
269            None => {
270                if self.sd.veh.veh_pt_type == vehicle::CONV {
271                    // If no EV / Hybrid components, no SOC considerations.
272                    (self.sd.veh.max_soc + self.sd.veh.min_soc) / 2.0
273                } else if self.sd.veh.veh_pt_type == vehicle::HEV {
274                    // ####################################
275                    // ### Charge Balancing Vehicle SOC ###
276                    // ####################################
277                    // Charge balancing SOC for HEV vehicle types. Iterating init_soc and comparing to final SOC.
278                    // Iterating until tolerance met or 30 attempts made.
279                    let mut init_soc = (self.sd.veh.max_soc + self.sd.veh.min_soc) / 2.0;
280                    let mut ess_2fuel_kwh = 1.0;
281                    while ess_2fuel_kwh > self.sd.veh.ess_to_fuel_ok_error
282                        && self.sd.hev_sim_count < self.sd.sim_params.sim_count_max
283                    {
284                        self.sd.hev_sim_count += 1;
285                        self.walk(init_soc, aux_in_kw_override.clone());
286                        let fuel_kj = (&self.sd.fs_kw_out_ach * self.sd.cyc.dt_s()).sum();
287                        let roadway_chg_kj =
288                            (&self.sd.roadway_chg_kw_out_ach * self.sd.cyc.dt_s()).sum();
289                        if (fuel_kj + roadway_chg_kj) > 0.0 {
290                            ess_2fuel_kwh = ((self.sd.soc[0] - self.sd.soc.last().unwrap())
291                                * self.sd.veh.ess_max_kwh
292                                * 3.6e3
293                                / (fuel_kj + roadway_chg_kj))
294                                .abs();
295                        } else {
296                            ess_2fuel_kwh = 0.0;
297                        }
298                        init_soc = min(1.0, max(0.0, *self.sd.soc.last().unwrap()));
299                    }
300                    init_soc
301                } else if self.sd.veh.veh_pt_type == vehicle::PHEV
302                    || self.sd.veh.veh_pt_type == vehicle::BEV
303                {
304                    // If EV, initializing initial SOC to maximum SOC.
305                    self.sd.veh.max_soc
306                } else {
307                    bail!("Failed to properly initialize SOC.");
308                }
309            }
310        };
311
312        self.walk(init_soc, aux_in_kw_override);
313
314        self.set_post_scalars()?;
315        Ok(())
316    }
317
318    pub fn walk(&mut self, init_soc: f64, aux_in_kw_override: Option<Array1<f64>>) {
319        self.init_for_step(init_soc, aux_in_kw_override);
320        while self.sd.i < self.sd.cyc.len() {
321            self.step().unwrap();
322        }
323    }
324
325    pub fn init_for_step(&mut self, init_soc: f64, aux_in_kw_override: Option<Array1<f64>>) {
326        self.history.push(self.state.clone()); // TODO: eventually make this dependent on `save_interval` usize per ALTRIOS
327        match &self.vehthrm.cabin_hvac_model {
328            CabinHvacModelTypes::Internal(hvac_mod) => {
329                self.hvac_model_history.push(hvac_mod.clone())
330            }
331            CabinHvacModelTypes::External => {}
332        }
333        self.sd.init_for_step(init_soc, aux_in_kw_override).unwrap();
334    }
335
336    pub fn set_speed_for_target_gap_using_idm(&mut self, i: usize) {
337        self.sd.set_speed_for_target_gap_using_idm(i);
338    }
339
340    pub fn set_speed_for_target_gap(&mut self, i: usize) {
341        self.sd.set_speed_for_target_gap(i);
342    }
343
344    pub fn step(&mut self) -> anyhow::Result<()> {
345        self.set_thermal_calcs(self.sd.i)?;
346        self.set_misc_calcs(self.sd.i);
347        self.set_comp_lims(self.sd.i)?;
348        self.set_power_calcs(self.sd.i)?;
349        self.set_ach_speed(self.sd.i)?;
350        self.set_hybrid_cont_calcs(self.sd.i)?;
351        self.set_fc_forced_state_rust(self.sd.i)?;
352        self.set_hybrid_cont_decisions(self.sd.i)?;
353        self.set_fc_power(self.sd.i)?;
354
355        self.sd.i += 1; // increment time step counter
356        self.history.push(self.state.clone());
357        match &self.vehthrm.cabin_hvac_model {
358            CabinHvacModelTypes::Internal(hvac_mod) => {
359                self.hvac_model_history.push(hvac_mod.clone());
360            }
361            CabinHvacModelTypes::External => {}
362        }
363        Ok(())
364    }
365
366    pub fn solve_step(&mut self, i: usize) -> anyhow::Result<()> {
367        self.sd.solve_step(i)
368    }
369
370    pub fn set_thermal_calcs(&mut self, i: usize) -> anyhow::Result<()> {
371        // most of the thermal equations are at [i-1] because the various thermally
372        // sensitive component efficiencies dependent on the [i] temperatures, but
373        // these are in turn dependent on [i-1] heat transfer processes
374        // verify that valid option is specified
375
376        if let Some(amb_te_deg_c) = &self.amb_te_deg_c {
377            self.state.amb_te_deg_c = amb_te_deg_c[i];
378        }
379
380        if let FcModelTypes::Internal(..) = &self.vehthrm.fc_model {
381            self.set_fc_thermal_calcs(i)?;
382        }
383
384        if let CabinHvacModelTypes::Internal(_) = &self.vehthrm.cabin_hvac_model {
385            self.set_cab_thermal_calcs(i)?;
386        }
387
388        if self.vehthrm.exhport_model == ComponentModelTypes::Internal {
389            self.set_exhport_thermal_calcs(i)?
390        }
391
392        if self.vehthrm.cat_model == ComponentModelTypes::Internal {
393            self.set_cat_thermal_calcs(i)?
394        }
395
396        if self.vehthrm.fc_model != FcModelTypes::External {
397            // Energy balance for fuel converter
398            self.state.fc_te_deg_c += (self.state.fc_qdot_kw
399                - self.state.fc_qdot_to_amb_kw
400                - self.state.fc_qdot_to_htr_kw)
401                / self.vehthrm.fc_c_kj__k
402                * self.sd.cyc.dt_s_at_i(i)
403        }
404        Ok(())
405    }
406
407    /// Solve fuel converter thermal behavior assuming convection parameters of sphere.
408    pub fn set_fc_thermal_calcs(&mut self, i: usize) -> anyhow::Result<()> {
409        // Constitutive equations for fuel converter
410        // calculation of adiabatic flame temperature
411        self.state.fc_te_adiabatic_deg_c = self.air.get_te_from_h(
412            ((1.0 + self.state.fc_lambda * self.sd.props.fuel_afr_stoich)
413                * self.air.get_h(self.state.amb_te_deg_c)?
414                + self.sd.props.get_fuel_lhv_kj_per_kg() * 1e3 * self.state.fc_lambda.min(1.0))
415                / (1.0 + self.state.fc_lambda * self.sd.props.fuel_afr_stoich),
416        )?;
417
418        // limited between 0 and 1, but should really not get near 1
419        self.state.fc_qdot_per_net_heat = (self.vehthrm.fc_coeff_from_comb
420            * (self.state.fc_te_adiabatic_deg_c - self.state.fc_te_deg_c))
421            .clamp(0.0, 1.0);
422
423        // heat generation
424        self.state.fc_qdot_kw = self.state.fc_qdot_per_net_heat
425            * (self.sd.fc_kw_in_ach[i - 1] - self.sd.fc_kw_out_ach[i - 1]);
426
427        // film temperature for external convection calculations
428        let fc_air_film_te_deg_c = 0.5 * (self.state.fc_te_deg_c + self.state.amb_te_deg_c);
429
430        // density * speed * diameter / dynamic viscosity
431        let fc_air_film_re = self.air.get_rho(fc_air_film_te_deg_c, None)
432            * self.sd.mps_ach[i - 1]
433            * self.vehthrm.fc_l
434            / self.air.get_mu(fc_air_film_te_deg_c)?;
435        // calculate heat transfer coeff. from engine to ambient [W / (m ** 2 * K)]
436        if self.sd.mps_ach[i - 1] < 1.0 {
437            // if stopped, scale based on thermostat opening and constant convection
438            self.state.fc_htc_to_amb = interpolate(
439                &self.state.fc_te_deg_c,
440                &Array1::from_vec(vec![
441                    self.vehthrm.tstat_te_sto_deg_c,
442                    self.vehthrm.tstat_te_fo_deg_c(),
443                ]),
444                &Array1::from_vec(vec![
445                    self.vehthrm.fc_htc_to_amb_stop,
446                    self.vehthrm.fc_htc_to_amb_stop * self.vehthrm.rad_eps,
447                ]),
448                false,
449            )
450            .with_context(|| format_dbg!())?
451        } else {
452            // Calculate heat transfer coefficient for sphere,
453            // from Incropera's Intro to Heat Transfer, 5th Ed., eq. 7.44
454            let fc_sphere_conv_params = get_sphere_conv_params(fc_air_film_re);
455            let fc_htc_to_amb_sphere = (fc_sphere_conv_params.0
456                * fc_air_film_re.powf(fc_sphere_conv_params.1))
457                * self.air.get_pr(fc_air_film_te_deg_c)?.powf(1.0 / 3.0)
458                * self.air.get_k(fc_air_film_te_deg_c)?
459                / self.vehthrm.fc_l;
460            self.state.fc_htc_to_amb = interpolate(
461                &self.state.fc_te_deg_c,
462                &Array1::from_vec(vec![
463                    self.vehthrm.tstat_te_sto_deg_c,
464                    self.vehthrm.tstat_te_fo_deg_c(),
465                ]),
466                &Array1::from_vec(vec![
467                    fc_htc_to_amb_sphere,
468                    fc_htc_to_amb_sphere * self.vehthrm.rad_eps,
469                ]),
470                false,
471            )
472            .with_context(|| format_dbg!())?
473        }
474
475        self.state.fc_qdot_to_amb_kw = self.state.fc_htc_to_amb
476            * 1e-3
477            * self.vehthrm.fc_area_ext()
478            * (self.state.fc_te_deg_c - self.state.amb_te_deg_c);
479        Ok(())
480    }
481
482    /// Solve cabin thermal behavior.
483    pub fn set_cab_thermal_calcs(&mut self, i: usize) -> anyhow::Result<()> {
484        if let CabinHvacModelTypes::Internal(hvac_model) = &mut self.vehthrm.cabin_hvac_model {
485            // flat plate model for isothermal, mixed-flow from Incropera and deWitt, Fundamentals of Heat and Mass
486            // Transfer, 7th Edition
487            let cab_te_film_ext_deg_c = 0.5 * (self.state.cab_te_deg_c + self.state.amb_te_deg_c);
488            let re_l = self.air.get_rho(cab_te_film_ext_deg_c, None)
489                * self.sd.mps_ach[i - 1]
490                * self.vehthrm.cab_l_length
491                / self.air.get_mu(cab_te_film_ext_deg_c)?;
492            let re_l_crit = 5.0e5; // critical Re for transition to turbulence
493
494            let nu_l_bar = if re_l < re_l_crit {
495                // equation 7.30
496                0.664 * re_l.powf(0.5) * self.air.get_pr(cab_te_film_ext_deg_c)?.powf(1.0 / 3.0)
497            } else {
498                // equation 7.38
499                let a = 871.0; // equation 7.39
500                (0.037 * re_l.powf(0.8) - a) * self.air.get_pr(cab_te_film_ext_deg_c)?
501            };
502
503            self.state.cab_qdot_to_amb_kw = if self.sd.mph_ach[i - 1] > 2.0 {
504                1e-3 * (self.vehthrm.cab_l_length * self.vehthrm.cab_l_width)
505                    / (1.0
506                        / (nu_l_bar * self.air.get_k(cab_te_film_ext_deg_c)?
507                            / self.vehthrm.cab_l_length)
508                        + self.vehthrm.cab_r_to_amb)
509                    * (self.state.cab_te_deg_c - self.state.amb_te_deg_c)
510            } else {
511                1e-3 * (self.vehthrm.cab_l_length * self.vehthrm.cab_l_width)
512                    / (1.0 / self.vehthrm.cab_htc_to_amb_stop + self.vehthrm.cab_r_to_amb)
513                    * (self.state.cab_te_deg_c - self.state.amb_te_deg_c)
514            };
515
516            let te_delta_vs_set_deg_c = self.state.cab_te_deg_c - hvac_model.te_set_deg_c;
517            let te_delta_vs_amb_deg_c = self.state.cab_te_deg_c - self.state.amb_te_deg_c;
518
519            if self.state.cab_te_deg_c <= hvac_model.te_set_deg_c + hvac_model.te_deadband_deg_c
520                && self.state.cab_te_deg_c >= hvac_model.te_set_deg_c - hvac_model.te_deadband_deg_c
521            {
522                // inside deadband; no hvac power is needed
523
524                self.state.cab_qdot_from_hvac_kw = 0.0;
525                hvac_model.i_cntrl_kw = 0.0; // reset to 0.0
526            } else {
527                hvac_model.p_cntrl_kw = hvac_model.p_cntrl_kw_per_deg_c * te_delta_vs_set_deg_c;
528                // integral control effort increases in magnitude by
529                // 1 time step worth of error
530                hvac_model.i_cntrl_kw += hvac_model.i_cntrl_kw_per_deg_c_scnds
531                    * te_delta_vs_set_deg_c
532                    * self.sd.cyc.dt_s_at_i(i);
533
534                hvac_model.d_cntrl_kw = hvac_model.d_cntrl_kj_per_deg_c
535                    * ((self.state.cab_te_deg_c - self.state.cab_prev_te_deg_c)
536                        / self.sd.cyc.dt_s_at_i(i));
537
538                // https://en.wikipedia.org/wiki/Coefficient_of_performance#Theoretical_performance_limits
539                // cop_ideal is t_h / (t_h - t_c) for heating
540                // cop_ideal is t_c / (t_h - t_c) for cooling
541
542                // divide-by-zero protection and realistic limit on COP
543                let cop_ideal = if te_delta_vs_amb_deg_c.abs() < 5.0 {
544                    // cabin is cooler than ambient + threshold
545                    (self.state.cab_te_deg_c + 273.15) / 5.0
546                } else {
547                    (self.state.cab_te_deg_c + 273.15) / te_delta_vs_amb_deg_c.abs()
548                };
549                hvac_model.cop = cop_ideal * hvac_model.frac_of_ideal_cop;
550                assert!(hvac_model.cop > 0.0);
551
552                if self.state.cab_te_deg_c > hvac_model.te_set_deg_c + hvac_model.te_deadband_deg_c
553                {
554                    // COOLING MODE; cabin is hotter than set point
555
556                    if hvac_model.i_cntrl_kw < 0.0 {
557                        // reset to switch from heating to cooling
558                        hvac_model.i_cntrl_kw = 0.0;
559                    }
560                    hvac_model.i_cntrl_kw = hvac_model.i_cntrl_kw.min(hvac_model.cntrl_max_kw);
561                    self.state.cab_qdot_from_hvac_kw =
562                        (-hvac_model.p_cntrl_kw - hvac_model.i_cntrl_kw - hvac_model.d_cntrl_kw)
563                            .max(-hvac_model.cntrl_max_kw);
564
565                    self.state.cab_hvac_pwr_aux_kw = (-self.state.cab_qdot_from_hvac_kw
566                        / hvac_model.cop)
567                        .min(hvac_model.pwr_max_aux_load_for_cooling_kw)
568                        .max(0.0);
569                    // correct if limit is exceeded
570                    self.state.cab_qdot_from_hvac_kw =
571                        -self.state.cab_hvac_pwr_aux_kw * hvac_model.cop;
572                } else {
573                    // HEATING MODE; cabin is colder than set point
574
575                    if hvac_model.i_cntrl_kw > 0.0 {
576                        // reset to switch from cooling to heating
577                        hvac_model.i_cntrl_kw = 0.0;
578                    }
579                    hvac_model.i_cntrl_kw = hvac_model.i_cntrl_kw.max(-hvac_model.cntrl_max_kw);
580
581                    self.state.cab_qdot_from_hvac_kw =
582                        (-hvac_model.p_cntrl_kw - hvac_model.i_cntrl_kw - hvac_model.d_cntrl_kw)
583                            .min(hvac_model.cntrl_max_kw);
584
585                    if hvac_model.use_fc_waste_heat {
586                        // limit heat transfer to be substantially less than what is physically possible
587                        // i.e. the engine can't drop below cabin temperature to heat the cabin
588                        self.state.cab_qdot_from_hvac_kw = self
589                            .state
590                            .cab_qdot_from_hvac_kw
591                            .min(
592                                (self.state.fc_te_deg_c - self.state.cab_te_deg_c)
593                                * 0.1 // so that it's substantially less
594                                * self.vehthrm.cab_c_kj__k
595                                    / self.sd.cyc.dt_s_at_i(i),
596                            )
597                            .max(0.0);
598                        self.state.fc_qdot_to_htr_kw = self.state.cab_qdot_from_hvac_kw;
599                        // TODO: think about what to do for PHEV, which needs careful consideration here
600                        // HEV probably also needs careful consideration
601                        // There needs to be an engine temperature (e.g. 60°C) below which the engine is forced on
602                        assert!(self.sd.veh.veh_pt_type != "BEV");
603                        // assume blower has negligible impact on aux load, may want to revise later
604                    } else {
605                        self.state.cab_hvac_pwr_aux_kw = (self.state.cab_qdot_from_hvac_kw
606                            / hvac_model.cop)
607                            .min(hvac_model.pwr_max_aux_load_for_cooling_kw)
608                            .max(0.0);
609                        self.state.cab_qdot_from_hvac_kw =
610                            self.state.cab_hvac_pwr_aux_kw * hvac_model.cop;
611                    }
612                }
613            }
614
615            self.state.cab_prev_te_deg_c = self.state.cab_te_deg_c;
616            self.state.cab_te_deg_c += (self.state.cab_qdot_from_hvac_kw
617                - self.state.cab_qdot_to_amb_kw)
618                / self.vehthrm.cab_c_kj__k
619                * self.sd.cyc.dt_s_at_i(i);
620        }
621        Ok(())
622    }
623
624    /// Solve exhport thermal behavior.
625    pub fn set_exhport_thermal_calcs(&mut self, i: usize) -> anyhow::Result<()> {
626        // lambda index may need adjustment, depending on how this ends up being modeled.
627        self.state.exh_mdot = self.sd.fs_kw_out_ach[i - 1] / self.sd.props.get_fuel_lhv_kj_per_kg()
628            * (1.0 + self.sd.props.fuel_afr_stoich * self.state.fc_lambda);
629        self.state.exh_hdot_kw = (1.0 - self.state.fc_qdot_per_net_heat)
630            * (self.sd.fc_kw_in_ach[i - 1] - self.sd.fc_kw_out_ach[i - 1]);
631
632        if self.state.exh_mdot > 5e-4 {
633            self.state.exhport_exh_te_in_deg_c = min(
634                self.air
635                    .get_te_from_h(self.state.exh_hdot_kw * 1e3 / self.state.exh_mdot)?,
636                self.state.fc_te_adiabatic_deg_c,
637            );
638            // when flow is small, assume inlet temperature is temporally constant
639            // so previous value is not overwritten
640        }
641
642        // calculate heat transfer coeff. from exhaust port to ambient [W / (m ** 2 * K)]
643        if (self.state.exhport_te_deg_c - self.state.fc_te_deg_c) > 0.0 {
644            // if exhaust port is hotter than ambient, make sure heat transfer cannot violate the second law
645            self.state.exhport_qdot_to_amb = min(
646                // nominal heat transfer to amb
647                self.vehthrm.exhport_ha_to_amb
648                    * (self.state.exhport_te_deg_c - self.state.fc_te_deg_c),
649                // max possible heat transfer to amb
650                self.vehthrm.exhport_c_kj__k
651                    * 1e3
652                    * (self.state.exhport_te_deg_c - self.state.fc_te_deg_c)
653                    / self.sd.cyc.dt_s_at_i(i),
654            );
655        } else {
656            // exhaust port cooler than the ambient
657            self.state.exhport_qdot_to_amb = max(
658                // nominal heat transfer to amb
659                self.vehthrm.exhport_ha_to_amb
660                    * (self.state.exhport_te_deg_c - self.state.fc_te_deg_c),
661                // max possible heat transfer to amb
662                self.vehthrm.exhport_c_kj__k
663                    * 1e3
664                    * (self.state.exhport_te_deg_c - self.state.fc_te_deg_c)
665                    / self.sd.cyc.dt_s_at_i(i),
666            );
667        }
668
669        if (self.state.exhport_exh_te_in_deg_c - self.state.exhport_te_deg_c) > 0.0 {
670            // exhaust hotter than exhaust port
671            self.state.exhport_qdot_from_exh = arrmin(&[
672                // nominal heat transfer to exhaust port
673                self.vehthrm.exhport_ha_int
674                    * (self.state.exhport_exh_te_in_deg_c - self.state.exhport_te_deg_c),
675                // max possible heat transfer from exhaust
676                self.state.exh_mdot
677                    * (self.air.get_h(self.state.exhport_exh_te_in_deg_c)?
678                        - self.air.get_h(self.state.exhport_te_deg_c)?),
679                // max possible heat transfer to exhaust port
680                self.vehthrm.exhport_c_kj__k
681                    * 1e3
682                    * (self.state.exhport_exh_te_in_deg_c - self.state.exhport_te_deg_c)
683                    / self.sd.cyc.dt_s_at_i(i),
684            ]);
685        } else {
686            // exhaust cooler than exhaust port
687            self.state.exhport_qdot_from_exh = arrmax(&[
688                // nominal heat transfer to exhaust port
689                self.vehthrm.exhport_ha_int
690                    * (self.state.exhport_exh_te_in_deg_c - self.state.exhport_te_deg_c),
691                // max possible heat transfer from exhaust
692                self.state.exh_mdot
693                    * (self.air.get_h(self.state.exhport_exh_te_in_deg_c)?
694                        - self.air.get_h(self.state.exhport_te_deg_c)?),
695                // max possible heat transfer to exhaust port
696                self.vehthrm.exhport_c_kj__k
697                    * 1e3
698                    * (self.state.exhport_exh_te_in_deg_c - self.state.exhport_te_deg_c)
699                    / self.sd.cyc.dt_s_at_i(i),
700            ]);
701        }
702
703        self.state.exhport_qdot_net =
704            self.state.exhport_qdot_from_exh - self.state.exhport_qdot_to_amb;
705        self.state.exhport_te_deg_c += self.state.exhport_qdot_net
706            / (self.vehthrm.exhport_c_kj__k * 1e3)
707            * self.sd.cyc.dt_s_at_i(i);
708        Ok(())
709    }
710
711    pub fn thermal_soak_walk(&mut self) -> anyhow::Result<()> {
712        self.sd.i = 1;
713        loop {
714            if self.sd.i < self.sd.cyc.len() {
715                break Ok(());
716            }
717            self.set_thermal_calcs(self.sd.i)?;
718            self.sd.i += 1;
719        }
720    }
721
722    /// Solve catalyst thermal behavior.
723    pub fn set_cat_thermal_calcs(&mut self, i: usize) -> anyhow::Result<()> {
724        // external or internal model handling catalyst thermal behavior
725
726        // Constitutive equations for catalyst
727        // catalyst film temperature for property calculation
728        let cat_te_ext_film_deg_c = 0.5 * (self.state.cat_te_deg_c + self.state.amb_te_deg_c);
729        // density * speed * diameter / dynamic viscosity
730        self.state.cat_re_ext = self.air.get_rho(cat_te_ext_film_deg_c, None)
731            * self.sd.mps_ach[i - 1]
732            * self.vehthrm.cat_l
733            / self.air.get_mu(cat_te_ext_film_deg_c)?;
734
735        // calculate heat transfer coeff. from cat to ambient [W / (m ** 2 * K)]
736        if self.sd.mps_ach[i - 1] < 1.0 {
737            // if stopped, scale based on constant convection
738            self.state.cat_htc_to_amb = self.vehthrm.cat_htc_to_amb_stop;
739        } else {
740            // if moving, scale based on speed dependent convection and thermostat opening
741            // Nusselt number coefficients from Incropera's Intro to Heat Transfer, 5th Ed., eq. 7.44
742            let cat_sphere_conv_params = get_sphere_conv_params(self.state.cat_re_ext);
743            self.state.fc_htc_to_amb = (cat_sphere_conv_params.0
744                * self.state.cat_re_ext.powf(cat_sphere_conv_params.1))
745                * self.air.get_pr(cat_te_ext_film_deg_c)?.powf(1.0 / 3.0)
746                * self.air.get_k(cat_te_ext_film_deg_c)?
747                / self.vehthrm.cat_l;
748        }
749
750        if (self.state.cat_te_deg_c - self.state.amb_te_deg_c) > 0.0 {
751            // cat hotter than ambient
752            self.state.cat_qdot_to_amb = min(
753                // nominal heat transfer to ambient
754                self.state.cat_htc_to_amb
755                    * self.vehthrm.cat_area_ext()
756                    * (self.state.cat_te_deg_c - self.state.amb_te_deg_c),
757                // max possible heat transfer to ambient
758                self.vehthrm.cat_c_kj__K
759                    * 1e3
760                    * (self.state.cat_te_deg_c - self.state.amb_te_deg_c)
761                    / self.sd.cyc.dt_s_at_i(i),
762            );
763        } else {
764            // ambient hotter than cat (less common)
765            self.state.cat_qdot_to_amb = max(
766                // nominal heat transfer to ambient
767                self.state.cat_htc_to_amb
768                    * self.vehthrm.cat_area_ext()
769                    * (self.state.cat_te_deg_c - self.state.amb_te_deg_c),
770                // max possible heat transfer to ambient
771                self.vehthrm.cat_c_kj__K
772                    * 1e3
773                    * (self.state.cat_te_deg_c - self.state.amb_te_deg_c)
774                    / self.sd.cyc.dt_s_at_i(i),
775            );
776        }
777
778        if self.state.exh_mdot > 5e-4 {
779            self.state.cat_exh_te_in_deg_c = min(
780                self.air.get_te_from_h(
781                    (self.state.exh_hdot_kw * 1e3 - self.state.exhport_qdot_from_exh)
782                        / self.state.exh_mdot,
783                )?,
784                self.state.fc_te_adiabatic_deg_c,
785            );
786            // when flow is small, assume inlet temperature is temporally constant
787            // so previous value is not overwritten
788        }
789
790        if (self.state.cat_exh_te_in_deg_c - self.state.cat_te_deg_c) > 0.0 {
791            // exhaust hotter than cat
792            self.state.cat_qdot_from_exh = min(
793                // limited by exhaust heat capacitance flow
794                self.state.exh_mdot
795                    * (self.air.get_h(self.state.cat_exh_te_in_deg_c)?
796                        - self.air.get_h(self.state.cat_te_deg_c)?),
797                // limited by catalyst thermal mass temperature change
798                self.vehthrm.cab_c_kj__k
799                    * 1e3
800                    * (self.state.cat_exh_te_in_deg_c - self.state.cat_te_deg_c)
801                    / self.sd.cyc.dt_s_at_i(i),
802            );
803        } else {
804            // cat hotter than exhaust (less common)
805            self.state.cat_qdot_from_exh = max(
806                // limited by exhaust heat capacitance flow
807                self.state.exh_mdot
808                    * (self.air.get_h(self.state.cat_exh_te_in_deg_c)?
809                        - self.air.get_h(self.state.cat_te_deg_c)?),
810                // limited by catalyst thermal mass temperature change
811                self.vehthrm.cat_c_kj__K
812                    * 1e3
813                    * (self.state.cat_exh_te_in_deg_c - self.state.cat_te_deg_c)
814                    / self.sd.cyc.dt_s_at_i(i),
815            );
816        }
817
818        // catalyst heat generation
819        self.state.cat_qdot = 0.0; // TODO: put something substantive here eventually
820
821        // net heat generetion/transfer in cat
822        self.state.cat_qdot_net =
823            self.state.cat_qdot + self.state.cat_qdot_from_exh - self.state.cat_qdot_to_amb;
824
825        self.state.cat_te_deg_c +=
826            self.state.cat_qdot_net * 1e-3 / self.vehthrm.cat_c_kj__K * self.sd.cyc.dt_s_at_i(i);
827        Ok(())
828    }
829
830    pub fn set_misc_calcs(&mut self, i: usize) {
831        // if cycle iteration is used, auxInKw must be re-zeroed to trigger the below if statement
832        // TODO: this is probably computationally expensive and was probably a workaround for numba
833        // figure out a way to not need this
834        if self.sd.aux_in_kw.slice(s![i..]).iter().all(|&x| x == 0.0) {
835            // if all elements after i-1 are zero, trigger default behavior; otherwise, use override value
836            if self.sd.veh.no_elec_aux {
837                self.sd.aux_in_kw[i] = self.sd.veh.aux_kw / self.sd.veh.alt_eff;
838            } else {
839                self.sd.aux_in_kw[i] = self.sd.veh.aux_kw;
840            }
841        }
842        self.sd.aux_in_kw[i] += self.state.cab_hvac_pwr_aux_kw;
843        // Is SOC below min threshold?
844        self.sd.reached_buff[i] =
845            self.sd.soc[i - 1] >= (self.sd.veh.min_soc + self.sd.veh.perc_high_acc_buf);
846
847        // Does the engine need to be on for low SOC or high acceleration
848        self.sd.high_acc_fc_on_tag[i] = self.sd.soc[i - 1] < self.sd.veh.min_soc
849            || (self.sd.high_acc_fc_on_tag[i - 1] && !(self.sd.reached_buff[i]));
850        self.sd.max_trac_mps[i] =
851            self.sd.mps_ach[i - 1] + (self.sd.veh.max_trac_mps2 * self.sd.cyc.dt_s_at_i(i));
852    }
853
854    pub fn set_comp_lims(&mut self, i: usize) -> anyhow::Result<()> {
855        self.sd.set_comp_lims(i)
856    }
857
858    pub fn set_power_calcs(&mut self, i: usize) -> anyhow::Result<()> {
859        self.sd.set_power_calcs(i)
860    }
861
862    pub fn set_ach_speed(&mut self, i: usize) -> anyhow::Result<()> {
863        self.sd.set_ach_speed(i)
864    }
865
866    pub fn set_hybrid_cont_calcs(&mut self, i: usize) -> anyhow::Result<()> {
867        self.sd.set_hybrid_cont_calcs(i)
868    }
869
870    pub fn set_fc_forced_state_rust(&mut self, i: usize) -> anyhow::Result<()> {
871        self.sd.set_fc_forced_state_rust(i)
872    }
873
874    pub fn set_hybrid_cont_decisions(&mut self, i: usize) -> anyhow::Result<()> {
875        self.sd.set_hybrid_cont_decisions(i)
876    }
877
878    pub fn set_fc_power(&mut self, i: usize) -> anyhow::Result<()> {
879        if self.sd.veh.fc_max_kw == 0.0 {
880            self.sd.fc_kw_out_ach[i] = 0.0;
881        } else if self.sd.veh.fc_eff_type == vehicle::H2FC {
882            self.sd.fc_kw_out_ach[i] = min(
883                self.sd.cur_max_fc_kw_out[i],
884                max(
885                    0.0,
886                    self.sd.mc_elec_kw_in_ach[i] + self.sd.aux_in_kw[i]
887                        - self.sd.ess_kw_out_ach[i]
888                        - self.sd.roadway_chg_kw_out_ach[i],
889                ),
890            );
891        } else if self.sd.veh.no_elec_sys
892            || self.sd.veh.no_elec_aux
893            || self.sd.high_acc_fc_on_tag[i]
894        {
895            self.sd.fc_kw_out_ach[i] = min(
896                self.sd.cur_max_fc_kw_out[i],
897                max(
898                    0.0,
899                    self.sd.trans_kw_in_ach[i] - self.sd.mc_mech_kw_out_ach[i]
900                        + self.sd.aux_in_kw[i],
901                ),
902            );
903        } else {
904            self.sd.fc_kw_out_ach[i] = min(
905                self.sd.cur_max_fc_kw_out[i],
906                max(
907                    0.0,
908                    self.sd.trans_kw_in_ach[i] - self.sd.mc_mech_kw_out_ach[i],
909                ),
910            );
911        }
912
913        if self.sd.veh.fc_max_kw == 0.0 {
914            self.sd.fc_kw_out_ach_pct[i] = 0.0;
915        } else {
916            self.sd.fc_kw_out_ach_pct[i] = self.sd.fc_kw_out_ach[i] / self.sd.veh.fc_max_kw;
917        }
918
919        if self.sd.fc_kw_out_ach[i] == 0.0 {
920            self.sd.fc_kw_in_ach[i] = 0.0;
921            self.sd.fc_kw_out_ach_pct[i] = 0.0;
922        } else {
923            if let FcModelTypes::Internal(fc_temp_eff_model, fc_temp_eff_comp) =
924                &self.vehthrm.fc_model
925            {
926                self.state.fc_eta_temp_coeff = match fc_temp_eff_model {
927                    FcTempEffModel::Linear(FcTempEffModelLinear {
928                        offset,
929                        slope,
930                        minimum,
931                    }) => max(*minimum, min(1.0, offset + slope * self.state.fc_te_deg_c)),
932                    FcTempEffModel::Exponential(FcTempEffModelExponential {
933                        offset,
934                        lag,
935                        minimum,
936                    }) => {
937                        match fc_temp_eff_comp {
938                            FcTempEffComponent::FuelConverter => (1.0
939                                - f64::exp(-1.0 / lag * (self.state.fc_te_deg_c - offset)))
940                            .max(*minimum),
941                            FcTempEffComponent::CatAndFC => {
942                                if self.state.cat_te_deg_c < self.vehthrm.cat_te_lightoff_deg_c {
943                                    let fc_eta_temp_coeff = (1.0
944                                        - f64::exp(-1.0 / lag * (self.state.fc_te_deg_c - offset)))
945                                    .max(*minimum);
946                                    // reduce efficiency to account for catalyst not being lit off
947                                    fc_eta_temp_coeff * self.vehthrm.cat_fc_eta_coeff
948                                } else {
949                                    1.0
950                                }
951                            }
952                            FcTempEffComponent::Catalyst => (1.0
953                                - f64::exp(-1.0 / lag * (self.state.cat_te_deg_c - offset)))
954                            .max(*minimum),
955                        }
956                    }
957                }
958            }
959
960            if self.sd.fc_kw_out_ach[i] == *self.sd.veh.input_kw_out_array.max()? {
961                self.sd.fc_kw_in_ach[i] = self.sd.fc_kw_out_ach[i]
962                    / (self.sd.veh.fc_eff_array.last().unwrap() * self.state.fc_eta_temp_coeff)
963            } else {
964                self.sd.fc_kw_in_ach[i] = self.sd.fc_kw_out_ach[i]
965                    / (self.sd.veh.fc_eff_array[max(
966                        1.0,
967                        (first_grtr(
968                            &self.sd.veh.fc_kw_out_array,
969                            min(
970                                self.sd.fc_kw_out_ach[i],
971                                self.sd.veh.input_kw_out_array.max()? - 0.001,
972                            ),
973                        )
974                        .unwrap()
975                            - 1) as f64,
976                    ) as usize])
977                    / self.state.fc_eta_temp_coeff
978            }
979        }
980
981        // fs out = fc in
982        self.sd.fs_kw_out_ach[i] = self.sd.fc_kw_in_ach[i];
983
984        self.sd.fs_kwh_out_ach[i] = self.sd.fs_kw_out_ach[i] * self.sd.cyc.dt_s_at_i(i) / 3.6e3;
985        Ok(())
986    }
987
988    pub fn set_time_dilation(&mut self, i: usize) -> anyhow::Result<()> {
989        self.sd.set_time_dilation(i)
990    }
991
992    pub fn set_post_scalars(&mut self) -> anyhow::Result<()> {
993        self.sd.set_post_scalars()
994    }
995}
996
997#[add_pyo3_api(
998    #[new]
999    #[pyo3(signature = (
1000        amb_te_deg_c=None,
1001        fc_te_deg_c_init=None,
1002        cab_te_deg_c_init=None,
1003        exhport_te_deg_c_init=None,
1004        cat_te_deg_c_init=None,
1005    ))]
1006    pub fn __new__(
1007        amb_te_deg_c: Option<f64>,
1008        fc_te_deg_c_init: Option<f64>,
1009        cab_te_deg_c_init: Option<f64>,
1010        exhport_te_deg_c_init: Option<f64>,
1011        cat_te_deg_c_init: Option<f64>,
1012    ) -> Self {
1013        Self::new(
1014            amb_te_deg_c,
1015            fc_te_deg_c_init,
1016            cab_te_deg_c_init,
1017            exhport_te_deg_c_init,
1018            cat_te_deg_c_init,
1019        )
1020    }
1021)]
1022#[allow(non_snake_case)]
1023#[derive(Deserialize, Serialize, Clone, Debug, PartialEq, HistoryVec)]
1024/// Struct containing thermal state variables for all thermal components
1025pub struct ThermalState {
1026    // fuel converter (engine) variables
1027    /// fuel converter (engine) temperature \[°C\]
1028    pub fc_te_deg_c: f64,
1029    /// fuel converter temperature efficiency correction
1030    pub fc_eta_temp_coeff: f64,
1031    /// fuel converter heat generation per total heat release minus shaft power
1032    pub fc_qdot_per_net_heat: f64,
1033    /// fuel converter heat generation \[kW\]
1034    pub fc_qdot_kw: f64,
1035    /// fuel converter convection to ambient \[kW\]
1036    pub fc_qdot_to_amb_kw: f64,
1037    /// fuel converter heat loss to heater core \[kW\]
1038    pub fc_qdot_to_htr_kw: f64,
1039    /// heat transfer coeff \[W / (m ** 2 * K)\] to amb after arbitration
1040    pub fc_htc_to_amb: f64,
1041    /// lambda (air/fuel ratio normalized w.r.t. stoich air/fuel ratio) -- 1 is reasonable default
1042    pub fc_lambda: f64,
1043    /// lambda-dependent adiabatic flame temperature
1044    pub fc_te_adiabatic_deg_c: f64,
1045
1046    // cabin (cab) variables
1047    /// cabin temperature \[°C\]
1048    pub cab_te_deg_c: f64,
1049    /// previous cabin temperature \[°C\]
1050    pub cab_prev_te_deg_c: f64,
1051    /// cabin solar load \[kw\]
1052    pub cab_qdot_solar_kw: f64,
1053    /// cabin convection to ambient \[kw\]
1054    pub cab_qdot_to_amb_kw: f64,
1055    /// heat transfer to cabin from hvac system
1056    pub cab_qdot_from_hvac_kw: f64,
1057    /// aux load from hvac
1058    pub cab_hvac_pwr_aux_kw: f64,
1059
1060    // exhaust variables
1061    /// exhaust mass flow rate \[kg/s\]
1062    pub exh_mdot: f64,
1063    /// exhaust enthalpy flow rate \[kw\]
1064    pub exh_hdot_kw: f64,
1065
1066    /// exhaust port (exhport) variables
1067    /// exhaust temperature at exhaust port inlet
1068    pub exhport_exh_te_in_deg_c: f64,
1069    /// heat transfer from exhport to amb \[kw\]
1070    pub exhport_qdot_to_amb: f64,
1071    /// catalyst temperature \[°C\]
1072    pub exhport_te_deg_c: f64,
1073    /// convection from exhaust to exhport \[W\]
1074    /// positive means exhport is receiving heat
1075    pub exhport_qdot_from_exh: f64,
1076    /// net heat generation in cat \[W\]
1077    pub exhport_qdot_net: f64,
1078
1079    // catalyst (cat) variables
1080    /// catalyst heat generation \[W\]
1081    pub cat_qdot: f64,
1082    /// catalytic converter convection coefficient to ambient \[W / (m ** 2 * K)\]
1083    pub cat_htc_to_amb: f64,
1084    /// heat transfer from catalyst to ambient \[W\]
1085    pub cat_qdot_to_amb: f64,
1086    /// catalyst temperature \[°C\]
1087    pub cat_te_deg_c: f64,
1088    /// exhaust temperature at cat inlet
1089    pub cat_exh_te_in_deg_c: f64,
1090    /// catalyst external reynolds number
1091    pub cat_re_ext: f64,
1092    /// convection from exhaust to cat \[W\]
1093    /// positive means cat is receiving heat
1094    pub cat_qdot_from_exh: f64,
1095    /// net heat generation in cat \[W\]
1096    pub cat_qdot_net: f64,
1097
1098    /// ambient temperature
1099    pub amb_te_deg_c: f64,
1100    #[serde(skip)]
1101    pub orphaned: bool,
1102}
1103
1104impl SerdeAPI for ThermalState {}
1105
1106impl ThermalState {
1107    pub fn new(
1108        amb_te_deg_c: Option<f64>,
1109        fc_te_deg_c_init: Option<f64>,
1110        cab_te_deg_c_init: Option<f64>,
1111        exhport_te_deg_c_init: Option<f64>,
1112        cat_te_deg_c_init: Option<f64>,
1113    ) -> Self {
1114        // Note default temperature is defined twice, see default()
1115        let default_te_deg_c = 22.0;
1116        let amb_te_deg_c = amb_te_deg_c.unwrap_or(default_te_deg_c);
1117        Self {
1118            amb_te_deg_c,
1119            fc_te_deg_c: fc_te_deg_c_init.unwrap_or(amb_te_deg_c),
1120            cab_te_deg_c: cab_te_deg_c_init.unwrap_or(amb_te_deg_c),
1121            cab_prev_te_deg_c: cab_te_deg_c_init.unwrap_or(amb_te_deg_c),
1122            exhport_te_deg_c: exhport_te_deg_c_init.unwrap_or(amb_te_deg_c),
1123            cat_te_deg_c: cat_te_deg_c_init.unwrap_or(amb_te_deg_c),
1124            // fc_te_adiabatic_deg_c // chad is pretty sure 'fc_te_adiabatic_deg_c' gets overridden in first time step
1125            ..Default::default()
1126        }
1127    }
1128}
1129
1130impl Default for ThermalState {
1131    fn default() -> Self {
1132        // Note default temperature is defined twice, see new()
1133        let default_te_deg_c = 22.0;
1134
1135        Self {
1136            fc_te_deg_c: default_te_deg_c, // overridden by new()
1137            fc_eta_temp_coeff: 0.0,
1138            fc_qdot_per_net_heat: 0.0,
1139            fc_qdot_kw: 0.0,
1140            fc_qdot_to_amb_kw: 0.0,
1141            fc_qdot_to_htr_kw: 0.0,
1142            fc_htc_to_amb: 0.0,
1143            fc_lambda: 1.0,
1144            fc_te_adiabatic_deg_c: default_te_deg_c, // this needs to be calculated, get Chad to revisit
1145
1146            cab_te_deg_c: default_te_deg_c, // overridden by new()
1147            cab_prev_te_deg_c: default_te_deg_c,
1148            cab_qdot_solar_kw: 0.0,
1149            cab_qdot_to_amb_kw: 0.0,
1150            cab_qdot_from_hvac_kw: 0.0,
1151            cab_hvac_pwr_aux_kw: 0.0,
1152
1153            exh_mdot: 0.0,
1154            exh_hdot_kw: 0.0,
1155
1156            exhport_exh_te_in_deg_c: default_te_deg_c,
1157            exhport_qdot_to_amb: 0.0,
1158            exhport_te_deg_c: default_te_deg_c, // overridden by new()
1159            exhport_qdot_from_exh: 0.0,
1160            exhport_qdot_net: 0.0,
1161
1162            cat_qdot: 0.0,
1163            cat_htc_to_amb: 0.0,
1164            cat_qdot_to_amb: 0.0,
1165            cat_te_deg_c: default_te_deg_c, // overridden by new()
1166            cat_exh_te_in_deg_c: default_te_deg_c,
1167            cat_re_ext: 0.0,
1168            cat_qdot_from_exh: 0.0,
1169            cat_qdot_net: 0.0,
1170            amb_te_deg_c: default_te_deg_c, // overridden by new()
1171
1172            orphaned: false,
1173        }
1174    }
1175}