1use 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 #[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 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 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 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 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 pub fn sim_drive_step(&mut self) -> anyhow::Result<()> {
91 self.step()
92 }
93 #[pyo3(name = "solve_step")]
94 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 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 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 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 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 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 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 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 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 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 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, 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 (self.sd.veh.max_soc + self.sd.veh.min_soc) / 2.0
273 } else if self.sd.veh.veh_pt_type == vehicle::HEV {
274 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 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()); 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; 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 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 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 pub fn set_fc_thermal_calcs(&mut self, i: usize) -> anyhow::Result<()> {
409 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 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 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 let fc_air_film_te_deg_c = 0.5 * (self.state.fc_te_deg_c + self.state.amb_te_deg_c);
429
430 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 if self.sd.mps_ach[i - 1] < 1.0 {
437 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 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 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 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; let nu_l_bar = if re_l < re_l_crit {
495 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 let a = 871.0; (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 self.state.cab_qdot_from_hvac_kw = 0.0;
525 hvac_model.i_cntrl_kw = 0.0; } else {
527 hvac_model.p_cntrl_kw = hvac_model.p_cntrl_kw_per_deg_c * te_delta_vs_set_deg_c;
528 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 let cop_ideal = if te_delta_vs_amb_deg_c.abs() < 5.0 {
544 (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 if hvac_model.i_cntrl_kw < 0.0 {
557 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 self.state.cab_qdot_from_hvac_kw =
571 -self.state.cab_hvac_pwr_aux_kw * hvac_model.cop;
572 } else {
573 if hvac_model.i_cntrl_kw > 0.0 {
576 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 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 * 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 assert!(self.sd.veh.veh_pt_type != "BEV");
603 } 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 pub fn set_exhport_thermal_calcs(&mut self, i: usize) -> anyhow::Result<()> {
626 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 }
641
642 if (self.state.exhport_te_deg_c - self.state.fc_te_deg_c) > 0.0 {
644 self.state.exhport_qdot_to_amb = min(
646 self.vehthrm.exhport_ha_to_amb
648 * (self.state.exhport_te_deg_c - self.state.fc_te_deg_c),
649 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 self.state.exhport_qdot_to_amb = max(
658 self.vehthrm.exhport_ha_to_amb
660 * (self.state.exhport_te_deg_c - self.state.fc_te_deg_c),
661 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 self.state.exhport_qdot_from_exh = arrmin(&[
672 self.vehthrm.exhport_ha_int
674 * (self.state.exhport_exh_te_in_deg_c - self.state.exhport_te_deg_c),
675 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 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 self.state.exhport_qdot_from_exh = arrmax(&[
688 self.vehthrm.exhport_ha_int
690 * (self.state.exhport_exh_te_in_deg_c - self.state.exhport_te_deg_c),
691 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 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 pub fn set_cat_thermal_calcs(&mut self, i: usize) -> anyhow::Result<()> {
724 let cat_te_ext_film_deg_c = 0.5 * (self.state.cat_te_deg_c + self.state.amb_te_deg_c);
729 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 if self.sd.mps_ach[i - 1] < 1.0 {
737 self.state.cat_htc_to_amb = self.vehthrm.cat_htc_to_amb_stop;
739 } else {
740 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 self.state.cat_qdot_to_amb = min(
753 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 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 self.state.cat_qdot_to_amb = max(
766 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 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 }
789
790 if (self.state.cat_exh_te_in_deg_c - self.state.cat_te_deg_c) > 0.0 {
791 self.state.cat_qdot_from_exh = min(
793 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 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 self.state.cat_qdot_from_exh = max(
806 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 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 self.state.cat_qdot = 0.0; 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 self.sd.aux_in_kw.slice(s![i..]).iter().all(|&x| x == 0.0) {
835 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 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 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 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 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)]
1024pub struct ThermalState {
1026 pub fc_te_deg_c: f64,
1029 pub fc_eta_temp_coeff: f64,
1031 pub fc_qdot_per_net_heat: f64,
1033 pub fc_qdot_kw: f64,
1035 pub fc_qdot_to_amb_kw: f64,
1037 pub fc_qdot_to_htr_kw: f64,
1039 pub fc_htc_to_amb: f64,
1041 pub fc_lambda: f64,
1043 pub fc_te_adiabatic_deg_c: f64,
1045
1046 pub cab_te_deg_c: f64,
1049 pub cab_prev_te_deg_c: f64,
1051 pub cab_qdot_solar_kw: f64,
1053 pub cab_qdot_to_amb_kw: f64,
1055 pub cab_qdot_from_hvac_kw: f64,
1057 pub cab_hvac_pwr_aux_kw: f64,
1059
1060 pub exh_mdot: f64,
1063 pub exh_hdot_kw: f64,
1065
1066 pub exhport_exh_te_in_deg_c: f64,
1069 pub exhport_qdot_to_amb: f64,
1071 pub exhport_te_deg_c: f64,
1073 pub exhport_qdot_from_exh: f64,
1076 pub exhport_qdot_net: f64,
1078
1079 pub cat_qdot: f64,
1082 pub cat_htc_to_amb: f64,
1084 pub cat_qdot_to_amb: f64,
1086 pub cat_te_deg_c: f64,
1088 pub cat_exh_te_in_deg_c: f64,
1090 pub cat_re_ext: f64,
1092 pub cat_qdot_from_exh: f64,
1095 pub cat_qdot_net: f64,
1097
1098 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 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 ..Default::default()
1126 }
1127 }
1128}
1129
1130impl Default for ThermalState {
1131 fn default() -> Self {
1132 let default_te_deg_c = 22.0;
1134
1135 Self {
1136 fc_te_deg_c: default_te_deg_c, 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, cab_te_deg_c: default_te_deg_c, 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, 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, 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, orphaned: false,
1173 }
1174 }
1175}