Entropyk/crates/solver/src/coupling.rs

436 lines
14 KiB
Rust
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

//! Thermal coupling between circuits for heat transfer.
//!
//! This module provides the infrastructure for modeling heat exchange between
//! independent fluid circuits. Thermal couplings represent heat exchangers
//! that transfer heat from a "hot" circuit to a "cold" circuit without
//! fluid mixing.
//!
//! ## Sign Convention
//!
//! Heat transfer Q > 0 means heat flows INTO the cold circuit (out of hot circuit).
//! This follows the convention that the cold circuit receives heat.
//!
//! ## Coupling Graph and Circular Dependencies
//!
//! Thermal couplings form a directed graph where:
//! - Nodes are circuits (CircuitId)
//! - Edges point from hot_circuit to cold_circuit (direction of heat flow)
//!
//! Circular dependencies occur when circuits mutually heat each other (A→B and B→A).
//! Circuits in circular dependencies must be solved simultaneously by the solver.
use entropyk_core::{Temperature, ThermalConductance};
use petgraph::algo::{is_cyclic_directed, kosaraju_scc};
use petgraph::graph::{DiGraph, NodeIndex};
use std::collections::HashMap;
use crate::system::CircuitId;
/// Thermal coupling between two circuits via a heat exchanger.
///
/// Heat flows from `hot_circuit` to `cold_circuit` proportional to the
/// temperature difference and thermal conductance (UA value).
#[derive(Debug, Clone, PartialEq)]
pub struct ThermalCoupling {
/// Circuit that supplies heat (higher temperature side).
pub hot_circuit: CircuitId,
/// Circuit that receives heat (lower temperature side).
pub cold_circuit: CircuitId,
/// Thermal conductance (UA) in W/K. Higher values = more heat transfer.
pub ua: ThermalConductance,
/// Efficiency factor (0.0 to 1.0). Default is 1.0 (no losses).
pub efficiency: f64,
}
impl ThermalCoupling {
/// Creates a new thermal coupling between two circuits.
///
/// # Arguments
///
/// * `hot_circuit` - Circuit at higher temperature (heat source)
/// * `cold_circuit` - Circuit at lower temperature (heat sink)
/// * `ua` - Thermal conductance in W/K
///
/// # Example
///
/// ```
/// use entropyk_solver::{ThermalCoupling, CircuitId};
/// use entropyk_core::ThermalConductance;
///
/// let coupling = ThermalCoupling::new(
/// CircuitId(0),
/// CircuitId(1),
/// ThermalConductance::from_watts_per_kelvin(1000.0),
/// );
/// ```
pub fn new(hot_circuit: CircuitId, cold_circuit: CircuitId, ua: ThermalConductance) -> Self {
Self {
hot_circuit,
cold_circuit,
ua,
efficiency: 1.0,
}
}
/// Sets the efficiency factor for the coupling.
///
/// Efficiency accounts for heat losses in the heat exchanger.
/// A value of 0.9 means 90% of theoretical heat is transferred.
pub fn with_efficiency(mut self, efficiency: f64) -> Self {
self.efficiency = efficiency.clamp(0.0, 1.0);
self
}
}
/// Computes heat transfer for a thermal coupling.
///
/// # Formula
///
/// Q = η × UA × (T_hot - T_cold)
///
/// Where:
/// - Q is the heat transfer rate (W), positive means heat INTO cold circuit
/// - η is the efficiency factor
/// - UA is the thermal conductance (W/K)
/// - T_hot, T_cold are temperatures (K)
///
/// # Sign Convention
///
/// - Q > 0: Heat flows from hot to cold (normal operation)
/// - Q = 0: No temperature difference
/// - Q < 0: Cold is hotter than hot (reverse flow, unusual)
///
/// # Example
///
/// ```
/// use entropyk_solver::{ThermalCoupling, CircuitId, compute_coupling_heat};
/// use entropyk_core::{Temperature, ThermalConductance};
///
/// let coupling = ThermalCoupling::new(
/// CircuitId(0),
/// CircuitId(1),
/// ThermalConductance::from_watts_per_kelvin(1000.0),
/// );
///
/// let t_hot = Temperature::from_kelvin(350.0);
/// let t_cold = Temperature::from_kelvin(300.0);
///
/// let q = compute_coupling_heat(&coupling, t_hot, t_cold);
/// assert!(q > 0.0, "Heat should flow from hot to cold");
/// ```
pub fn compute_coupling_heat(
coupling: &ThermalCoupling,
t_hot: Temperature,
t_cold: Temperature,
) -> f64 {
coupling.efficiency
* coupling.ua.to_watts_per_kelvin()
* (t_hot.to_kelvin() - t_cold.to_kelvin())
}
/// Builds a coupling graph for dependency analysis.
///
/// Returns a directed graph where:
/// - Nodes are CircuitIds present in any coupling
/// - Edges point from hot_circuit to cold_circuit
fn build_coupling_graph(couplings: &[ThermalCoupling]) -> DiGraph<CircuitId, ()> {
let mut graph = DiGraph::new();
let mut circuit_to_node: HashMap<CircuitId, NodeIndex> = HashMap::new();
for coupling in couplings {
// Add hot_circuit node if not present
let hot_node = *circuit_to_node
.entry(coupling.hot_circuit)
.or_insert_with(|| graph.add_node(coupling.hot_circuit));
// Add cold_circuit node if not present
let cold_node = *circuit_to_node
.entry(coupling.cold_circuit)
.or_insert_with(|| graph.add_node(coupling.cold_circuit));
// Add directed edge: hot -> cold
graph.add_edge(hot_node, cold_node, ());
}
graph
}
/// Checks if the coupling graph contains circular dependencies.
///
/// Circular dependencies occur when circuits are mutually thermally coupled
/// (e.g., A heats B, and B heats A). When circular dependencies exist,
/// the solver must solve those circuits simultaneously rather than sequentially.
///
/// # Example
///
/// ```
/// use entropyk_solver::{ThermalCoupling, CircuitId, has_circular_dependencies};
/// use entropyk_core::ThermalConductance;
///
/// // No circular dependency: A → B → C
/// let couplings = vec![
/// ThermalCoupling::new(CircuitId(0), CircuitId(1), ThermalConductance::from_watts_per_kelvin(100.0)),
/// ThermalCoupling::new(CircuitId(1), CircuitId(2), ThermalConductance::from_watts_per_kelvin(100.0)),
/// ];
/// assert!(!has_circular_dependencies(&couplings));
///
/// // Circular dependency: A → B and B → A
/// let couplings_circular = vec![
/// ThermalCoupling::new(CircuitId(0), CircuitId(1), ThermalConductance::from_watts_per_kelvin(100.0)),
/// ThermalCoupling::new(CircuitId(1), CircuitId(0), ThermalConductance::from_watts_per_kelvin(100.0)),
/// ];
/// assert!(has_circular_dependencies(&couplings_circular));
/// ```
pub fn has_circular_dependencies(couplings: &[ThermalCoupling]) -> bool {
if couplings.is_empty() {
return false;
}
let graph = build_coupling_graph(couplings);
is_cyclic_directed(&graph)
}
/// Returns groups of circuits that must be solved simultaneously.
///
/// Groups are computed using strongly connected components (SCC) analysis
/// of the coupling graph. Circuits in the same SCC have circular thermal
/// dependencies and must be solved together.
///
/// # Returns
///
/// A vector of vectors, where each inner vector contains CircuitIds that
/// must be solved simultaneously. Single-element vectors indicate circuits
/// that can be solved independently (in topological order).
///
/// # Example
///
/// ```
/// use entropyk_solver::{ThermalCoupling, CircuitId, coupling_groups};
/// use entropyk_core::ThermalConductance;
///
/// // A → B, B and C independent
/// let couplings = vec![
/// ThermalCoupling::new(CircuitId(0), CircuitId(1), ThermalConductance::from_watts_per_kelvin(100.0)),
/// ];
/// let groups = coupling_groups(&couplings);
/// // Groups will contain individual circuits since there's no cycle
/// ```
pub fn coupling_groups(couplings: &[ThermalCoupling]) -> Vec<Vec<CircuitId>> {
if couplings.is_empty() {
return Vec::new();
}
let graph = build_coupling_graph(couplings);
let sccs = kosaraju_scc(&graph);
sccs.into_iter()
.map(|node_indices| node_indices.into_iter().map(|idx| graph[idx]).collect())
.collect()
}
#[cfg(test)]
mod tests {
use super::*;
use approx::assert_relative_eq;
fn make_coupling(hot: u8, cold: u8, ua_w_per_k: f64) -> ThermalCoupling {
ThermalCoupling::new(
CircuitId(hot),
CircuitId(cold),
ThermalConductance::from_watts_per_kelvin(ua_w_per_k),
)
}
#[test]
fn test_thermal_coupling_creation() {
let coupling = ThermalCoupling::new(
CircuitId(0),
CircuitId(1),
ThermalConductance::from_watts_per_kelvin(1000.0),
);
assert_eq!(coupling.hot_circuit, CircuitId(0));
assert_eq!(coupling.cold_circuit, CircuitId(1));
assert_relative_eq!(coupling.ua.to_watts_per_kelvin(), 1000.0, epsilon = 1e-10);
assert_relative_eq!(coupling.efficiency, 1.0, epsilon = 1e-10);
}
#[test]
fn test_thermal_coupling_with_efficiency() {
let coupling = ThermalCoupling::new(
CircuitId(0),
CircuitId(1),
ThermalConductance::from_watts_per_kelvin(1000.0),
)
.with_efficiency(0.85);
assert_relative_eq!(coupling.efficiency, 0.85, epsilon = 1e-10);
}
#[test]
fn test_efficiency_clamped() {
let coupling = make_coupling(0, 1, 100.0).with_efficiency(1.5);
assert_relative_eq!(coupling.efficiency, 1.0, epsilon = 1e-10);
let coupling = make_coupling(0, 1, 100.0).with_efficiency(-0.5);
assert_relative_eq!(coupling.efficiency, 0.0, epsilon = 1e-10);
}
#[test]
fn test_compute_coupling_heat_positive() {
let coupling = make_coupling(0, 1, 1000.0);
let t_hot = Temperature::from_kelvin(350.0);
let t_cold = Temperature::from_kelvin(300.0);
let q = compute_coupling_heat(&coupling, t_hot, t_cold);
// Q = 1.0 * 1000 * (350 - 300) = 50000 W
assert_relative_eq!(q, 50000.0, epsilon = 1e-10);
assert!(q > 0.0, "Heat should be positive (into cold circuit)");
}
#[test]
fn test_compute_coupling_heat_zero() {
let coupling = make_coupling(0, 1, 1000.0);
let t_hot = Temperature::from_kelvin(300.0);
let t_cold = Temperature::from_kelvin(300.0);
let q = compute_coupling_heat(&coupling, t_hot, t_cold);
assert_relative_eq!(q, 0.0, epsilon = 1e-10);
}
#[test]
fn test_compute_coupling_heat_negative() {
let coupling = make_coupling(0, 1, 1000.0);
let t_hot = Temperature::from_kelvin(280.0);
let t_cold = Temperature::from_kelvin(300.0);
let q = compute_coupling_heat(&coupling, t_hot, t_cold);
// Q = 1000 * (280 - 300) = -20000 W (reverse flow)
assert_relative_eq!(q, -20000.0, epsilon = 1e-10);
assert!(q < 0.0, "Heat should be negative (reverse flow)");
}
#[test]
fn test_compute_coupling_heat_with_efficiency() {
let coupling = make_coupling(0, 1, 1000.0).with_efficiency(0.9);
let t_hot = Temperature::from_kelvin(350.0);
let t_cold = Temperature::from_kelvin(300.0);
let q = compute_coupling_heat(&coupling, t_hot, t_cold);
// Q = 0.9 * 1000 * 50 = 45000 W
assert_relative_eq!(q, 45000.0, epsilon = 1e-10);
}
#[test]
fn test_energy_conservation() {
// For two circuits coupled, Q_hot = -Q_cold
// This means the heat leaving hot circuit equals heat entering cold circuit
let coupling = make_coupling(0, 1, 1000.0);
let t_hot = Temperature::from_kelvin(350.0);
let t_cold = Temperature::from_kelvin(300.0);
let q_into_cold = compute_coupling_heat(&coupling, t_hot, t_cold);
let q_out_of_hot = -q_into_cold; // By convention
// Heat into cold = - (heat out of hot)
assert_relative_eq!(q_into_cold, -q_out_of_hot, epsilon = 1e-10);
assert!(q_into_cold > 0.0, "Cold circuit receives heat");
assert!(q_out_of_hot < 0.0, "Hot circuit loses heat");
}
#[test]
fn test_no_circular_dependency() {
// Linear chain: A → B → C
let couplings = vec![make_coupling(0, 1, 100.0), make_coupling(1, 2, 100.0)];
assert!(!has_circular_dependencies(&couplings));
}
#[test]
fn test_circular_dependency_detection() {
// Mutual: A → B and B → A
let couplings = vec![make_coupling(0, 1, 100.0), make_coupling(1, 0, 100.0)];
assert!(has_circular_dependencies(&couplings));
}
#[test]
fn test_circular_dependency_complex() {
// Triangle: A → B → C → A
let couplings = vec![
make_coupling(0, 1, 100.0),
make_coupling(1, 2, 100.0),
make_coupling(2, 0, 100.0),
];
assert!(has_circular_dependencies(&couplings));
}
#[test]
fn test_empty_couplings_no_cycle() {
let couplings: Vec<ThermalCoupling> = vec![];
assert!(!has_circular_dependencies(&couplings));
}
#[test]
fn test_single_coupling_no_cycle() {
let couplings = vec![make_coupling(0, 1, 100.0)];
assert!(!has_circular_dependencies(&couplings));
}
#[test]
fn test_coupling_groups_no_cycle() {
// A → B, C independent
let couplings = vec![make_coupling(0, 1, 100.0)];
let groups = coupling_groups(&couplings);
// With no cycles, each circuit is its own group
assert_eq!(groups.len(), 2);
// Each group should have exactly one circuit
for group in &groups {
assert_eq!(group.len(), 1);
}
// Collect all circuit IDs
let all_circuits: std::collections::HashSet<CircuitId> =
groups.iter().flat_map(|g| g.iter().copied()).collect();
assert!(all_circuits.contains(&CircuitId(0)));
assert!(all_circuits.contains(&CircuitId(1)));
}
#[test]
fn test_coupling_groups_with_cycle() {
// A ↔ B (mutual), C → D
let couplings = vec![
make_coupling(0, 1, 100.0),
make_coupling(1, 0, 100.0),
make_coupling(2, 3, 100.0),
];
let groups = coupling_groups(&couplings);
// Should have 3 groups: [A, B] as one, C as one, D as one
assert_eq!(groups.len(), 3);
// Find the group with 2 circuits (A and B)
let large_group: Vec<&Vec<CircuitId>> = groups.iter().filter(|g| g.len() == 2).collect();
assert_eq!(large_group.len(), 1);
let ab_group = large_group[0];
assert!(ab_group.contains(&CircuitId(0)));
assert!(ab_group.contains(&CircuitId(1)));
}
#[test]
fn test_coupling_groups_empty() {
let couplings: Vec<ThermalCoupling> = vec![];
let groups = coupling_groups(&couplings);
assert!(groups.is_empty());
}
}