Coverage for pyrc\core\resistors.py: 60%

188 statements  

« prev     ^ index     » next       coverage.py v7.14.1, created at 2026-06-29 15:57 +0200

1# ------------------------------------------------------------------------------- 

2# Copyright (C) 2026 Joel Kimmich, Tim Jourdan 

3# ------------------------------------------------------------------------------ 

4# License 

5# This file is part of PyRC, distributed under GPL-3.0-or-later. 

6# ------------------------------------------------------------------------------ 

7 

8from __future__ import annotations 

9 

10import warnings 

11from typing import TYPE_CHECKING, Any 

12 

13import numpy as np 

14from scipy.constants import Stefan_Boltzmann 

15from sympy import Expr 

16 

17from pyrc.core.components.capacitor import Capacitor 

18from pyrc.core.components.resistor import Resistor 

19from pyrc.core.components.templates import Fluid, Solid 

20from pyrc.core.heat_transfer import alpha_forced_convection_in_pipe 

21from pyrc.tools.functions import check_type, is_set, return_type 

22 

23if TYPE_CHECKING: 

24 from pyrc.core.components.templates import TemperatureNode 

25 from pyrc.core.inputs import BoundaryCondition 

26 from pyrc.core.nodes import ChannelNode, MassFlowNode, Node 

27 

28np.seterr(divide="ignore") 

29 

30 

31class MassTransport(Resistor): 

32 def __init__(self): 

33 """ 

34 Represents the resistance caused by mass transfer between two :class:`MassFlowNode` s. 

35 

36 The resistance is calculated automatically. 

37 

38 Be aware that this :class:`Resistor` doesn't care about Courant number at all. This has to be checked in the Handler 

39 that starts the simulation. 

40 """ 

41 super().__init__(resistance=np.nan) 

42 self.__source: MassFlowNode | Any = None 

43 self.__sink: MassFlowNode | Any = None 

44 self._volume_flow: float | Any = None 

45 

46 @property 

47 def source(self) -> MassFlowNode: 

48 if self.__source is None: 

49 raise AttributeError("Source node has not been set yet.") 

50 return self.__source 

51 

52 @source.setter 

53 def source(self, value: MassFlowNode): 

54 self.__source = value 

55 

56 @property 

57 def sink(self) -> MassFlowNode: 

58 if self.__sink is None: 

59 raise AttributeError("Sink node has not been set yet.") 

60 return self.__sink 

61 

62 @sink.setter 

63 def sink(self, value: MassFlowNode): 

64 self.__sink = value 

65 

66 @property 

67 def guess_volume_flow(self): 

68 return self._volume_flow 

69 

70 @property 

71 def volume_flow(self): 

72 from pyrc.core.nodes import MassFlowNode 

73 

74 if self._volume_flow is None: 

75 # create volume flow using one connected MassFlowNode 

76 for node in self.neighbours: 

77 if isinstance(node, MassFlowNode): 

78 node.propagate_flow() 

79 break 

80 return self._volume_flow 

81 

82 @property 

83 def resistance(self) -> float | int: 

84 from pyrc.core.nodes import FlowBoundaryCondition 

85 

86 if not is_set(self._resistance): 

87 # Get volume flow using source and sink of neighbour nodes 

88 nodes: list = self.nodes 

89 assert len(nodes) == 2 

90 if isinstance(nodes[0], FlowBoundaryCondition): 

91 mfn = nodes[1] 

92 else: 

93 mfn = nodes[0] 

94 resistance = 1 / np.float64(self.volume_flow * mfn.material.density * mfn.material.heat_capacity) 

95 

96 self._resistance = resistance 

97 return self._resistance 

98 

99 # def make_velocity_direction(self): 

100 # """ 

101 # Trigger the process in a `MassFlowNode` neighbour of ``self``. 

102 # """ 

103 # from pyrc.core.nodes import MassFlowNode 

104 # for neighbour in self.neighbours: 

105 # if isinstance(neighbour, MassFlowNode): 

106 # neighbour.make_velocity_direction() 

107 # return None # break out of the method 

108 # # raise error otherwise 

109 # raise ConnectionError(f"The {type(self)} isn't connected to any MassFlowNode. Damn!") 

110 

111 

112class CombinedResistor(Resistor): 

113 def __init__( 

114 self, 

115 resistance: float | int | np.number | Expr = np.nan, 

116 htc: float | int | np.number | Expr = np.nan, 

117 heat_conduction=True, 

118 heat_transfer=True, 

119 ): 

120 """ 

121 Automated version of `Resistor`, calculating its resistance based on connected objects. 

122 

123 The algorithms of this class rely on geometric representations and assume a thermal problem. For most of the 

124 algorithms the connected capacities and boundary conditions must inherit from 

125 :class:`~pyrc.core.components.templates.Cell` or at least :class:`Geometric`\\. 

126 Using the geometric information let the algorithm figure out which areas/lengths influence the heat transfer and 

127 thermal conduction. If you connect a :class:`~pyrc.core.inputs.BoundaryCondition` (not inheriting from 

128 Cell/Geometric) you still can use this class but must define the direction to this BoundaryCondition 

129 manually using ``manual_directions`` and ``set_direction()`` of :class:`TemperatureNode` ( 

130 Capacitor/BoundaryCondition). 

131 

132 Parameters 

133 ---------- 

134 resistance : float | int | np.number | sympy.Expr, default=np.nan 

135 The resistance of self. If set, no algorithm is used (and it would be preferable to use :class:`Resistor` class 

136 instead). 

137 htc : float | int | np.number | sympy.Expr, default=np.nan 

138 Heat transfer coefficient that is used, if no other HTC was found (in boundary conditions). 

139 If not set, an initial value of 5 is used (raising a warning). 

140 heat_conduction : bool, default=True 

141 Switch on/off heat conductivity. 

142 heat_transfer : bool, default=True 

143 Switch on/off heat transfer. 

144 

145 See Also 

146 -------- 

147 Resistor : The basic Resistor class without the automatics. 

148 """ 

149 super().__init__(resistance) 

150 self.heat_conduction: bool = heat_conduction 

151 self.heat_transfer: bool = heat_transfer 

152 self.__htc = htc 

153 

154 @property 

155 def htc(self): 

156 if not is_set(self.__htc): 

157 warnings.warn("Warning: Initial HTC value of 5 is used.") 

158 self.__htc = 5 

159 return self.__htc 

160 

161 @property 

162 def heat_transfer_coefficient(self): 

163 return self.htc 

164 

165 @property 

166 def resistance(self) -> float | int | Expr: 

167 """ 

168 Determines the resistance accordingly to the nodes the resistor is connected to. 

169 

170 To get this, look at the pictures of Joel Kimmich from 13.8.2025 

171 

172 Returns 

173 ------- 

174 

175 """ 

176 if is_set(self._resistance): 

177 return self._resistance 

178 resistance = np.float64(0) 

179 

180 from pyrc.core.inputs import BoundaryCondition, FluidBoundaryCondition, SolidBoundaryCondition 

181 from pyrc.core.nodes import ChannelNode, Node 

182 

183 if len(self.direct_connected_node_templates) == 2: 

184 if all(isinstance(neighbour, Node) for neighbour in self.neighbours): 

185 # no BoundaryCondition 

186 if ( 

187 all(isinstance(neighbour.material, Solid) for neighbour in self.neighbours) 

188 or all(isinstance(neighbour.material, Fluid) for neighbour in self.neighbours) 

189 and self.heat_conduction 

190 ): 

191 # heat conduction in both nodes (also for two ChannelNodes) 

192 node: Node 

193 for i, node in enumerate(self.neighbours): 

194 other_node = self.neighbours[(i + 1) % 2] 

195 resistance += node.get_conduction_length(other_node) / ( 

196 node.material.thermal_conductivity * node.get_area(other_node) 

197 ) 

198 elif check_type(self.neighbours, ChannelNode, Node): 

199 # ChannelNode to Solid: HeatTransfer AND HeatConduction 

200 channel_node, node = return_type(self.neighbours, ChannelNode) 

201 assert isinstance(node.material, Solid) 

202 if self.heat_transfer: 

203 # calculate heat transfer in pipe using Nusselt correlation(s) 

204 resistance += resistance_channel_heat_transfer(channel_node, node) 

205 if self.heat_conduction: 

206 resistance += node.get_conduction_length(channel_node) / ( 

207 node.material.thermal_conductivity * node.get_area(channel_node) 

208 ) 

209 else: 

210 # solid to fluid: heat transfer in fluid and conduction in solid 

211 solid, fluid = return_type(self.neighbours, Solid, [n.material for n in self.neighbours]) 

212 if self.heat_transfer: 

213 effective_area = min(solid.get_area(fluid), fluid.get_area(solid)) 

214 resistance += 1 / (self.htc * effective_area) 

215 if self.heat_conduction: 

216 resistance += solid.get_conduction_length(fluid) / ( 

217 solid.material.thermal_conductivity * solid.get_area(fluid) 

218 ) 

219 elif check_type(self.neighbours, BoundaryCondition, Node): 

220 # 1 BC and 1 Node 

221 if check_type(self.neighbours, FluidBoundaryCondition, Node): 

222 # 1 FluidBC and 1 Node 

223 bc, node = return_type(self.neighbours, FluidBoundaryCondition) 

224 if isinstance(node.material, Solid): 

225 if self.heat_transfer: 

226 resistance += resistance_bc_heat_transfer(bc, node) 

227 if self.heat_conduction: 

228 # for both: (both Fluid) and (Fluid BC and Solid Node) 

229 resistance += node.get_conduction_length(bc) / ( 

230 node.material.thermal_conductivity * node.get_area(bc) 

231 ) 

232 elif check_type(self.neighbours, SolidBoundaryCondition, Node): 

233 # 1 SolidBC and 1 Node 

234 node, bc = return_type(self.neighbours, Node) 

235 effective_area = node.get_area(bc) 

236 if isinstance(node.material, Fluid) and self.heat_transfer: 

237 resistance += 1 / (self.htc * effective_area) 

238 # no heat conduction if node is Fluid (is already included in heat transfer) 

239 elif self.heat_conduction: 

240 # like two solids, but only node is used (so conduction in BC is infinite) 

241 resistance += node.get_conduction_length(bc) / ( 

242 node.material.thermal_conductivity * node.get_area(bc) 

243 ) 

244 else: 

245 # BC is undefined -> resistance should be given 

246 assert check_type(self.neighbours, BoundaryCondition, Node), ( 

247 "You must use FluidBoundaryCondition or SolidBoundaryCondition or set a resistance/htc value." 

248 ) 

249 bc, node = return_type(self.neighbours, BoundaryCondition) 

250 # TODO: Maybe a BC shouldn't get an HTC value because it depends on the cell / the cell AND the BC. 

251 if not is_set(bc.htc): 

252 htc = self.htc 

253 warnings.warn("Resistor.htc was used instead of BoundaryCondition.htc.") 

254 else: 

255 htc = bc.htc 

256 warnings.warn("You might consider using FluidBoundaryCondition instead of BoundaryCondition.") 

257 resistance += 1 / (htc * node.get_area(bc)) 

258 elif all(isinstance(neighbour, BoundaryCondition) for neighbour in self.neighbours): 

259 raise ValueError("Two BoundaryConditions cannot be connected directly to each other.") 

260 else: 

261 raise ValueError( 

262 f"This combination isn't implemented: {self.neighbours}\n{[type(n) for n in self.neighbours]}" 

263 ) 

264 elif len(self.direct_connected_node_templates) == 0: 

265 # Only connected to other resistors: 

266 # the resistance must be given. But this was already checked, so raise an Error 

267 raise ValueError("If a Resistor is between two others, its resistance must be given!") 

268 else: 

269 # one Node, one/multiple Resistor/s -> only the resistance of this Resistor is returned, no equivalent 

270 # resistance! 

271 node: Capacitor = self.direct_connected_node_templates[0] 

272 other_node: Capacitor = self.get_connected_node(node) 

273 if isinstance(node, BoundaryCondition): 

274 assert isinstance(other_node, Node) 

275 if isinstance(node, FluidBoundaryCondition): 

276 if isinstance(other_node.material, Fluid): 

277 # FluidBC to Fluid Node -> Own resistance is zero, because BC don't have mass/volume 

278 resistance += np.float64(0) 

279 elif self.heat_transfer: 

280 # FluidBC to Solid Node -> HeatTransfer 

281 resistance += resistance_bc_heat_transfer(node, other_node) 

282 elif isinstance(node, SolidBoundaryCondition): 

283 if isinstance(other_node.material, Fluid) and self.heat_transfer: 

284 # SolidBC to Fluid -> HeatTransfer 

285 resistance += resistance_bc_heat_transfer(node, other_node) 

286 elif isinstance(node, ChannelNode) and self.heat_transfer: 

287 # ChannelNode should only be connected to a Solid Node over a HeatConduction Resistor 

288 assert isinstance(other_node, Node) and isinstance(other_node.material, Solid) 

289 # HeatTransfer in a round channel 

290 resistance += resistance_channel_heat_transfer(node, other_node) 

291 else: 

292 # one Node (no BC) and one/multiple resistors -> only heat mechanism of self is calculated 

293 node: Node 

294 if self.heat_conduction and ( 

295 isinstance(node.material, Solid) 

296 or ( 

297 isinstance(node.material, Fluid) 

298 and ( 

299 isinstance(other_node, FluidBoundaryCondition) 

300 or (isinstance(other_node, Node) and isinstance(other_node.material, Fluid)) 

301 ) 

302 ) 

303 ): 

304 resistance = node.get_conduction_length(other_node) / ( 

305 node.material.thermal_conductivity * node.get_area(other_node) 

306 ) 

307 elif self.heat_transfer: 

308 # self.material is Fluid and other_node is solid 

309 # Heat transfer 

310 if isinstance(other_node, SolidBoundaryCondition): 

311 resistance += resistance_bc_heat_transfer(other_node, node) 

312 else: 

313 resistance += 1 / (self.htc * node.get_area(other_node)) 

314 return resistance 

315 

316 

317class HeatConduction(CombinedResistor): 

318 def __init__(self, resistance=np.nan): 

319 """ 

320 Represents the resistance caused by heat conduction. 

321 

322 If the nodes, where the heat conduction takes place, differ in their material (Solid and Fluid) the heat 

323 conduction is set to 0 (the resistance is set to np.inf), because the heat conduction is included in 

324 HeatTransfer. So calculating it also in HeatConduction it would be taken into account twice. 

325 So: Do not forget to create a `HeatTransfer` `Resistor` between such nodes! 

326 

327 Parameters 

328 ---------- 

329 resistance : float, default=np.nan 

330 The resistance. If set, it will not be calculated. 

331 """ 

332 super().__init__(resistance, heat_conduction=True, heat_transfer=False) 

333 

334 @property 

335 def htc(self): 

336 return np.nan 

337 

338 

339class HeatTransfer(CombinedResistor): 

340 def __init__(self, resistance=np.nan, htc=np.nan): 

341 """ 

342 Represents the resistance caused by heat transfer between a solid and a fluid. 

343 

344 Parameters 

345 ---------- 

346 resistance : float, default=np.nan 

347 The resistance. If set, it will not be calculated. 

348 """ 

349 super().__init__(resistance, htc=htc, heat_transfer=True, heat_conduction=False) 

350 

351 # def single_resistance(self, node: FluidBoundaryCondition | Node): 

352 # """ 

353 # Returns the resistance using the passed node. 

354 # 

355 # The passed node must have a Fluid as material or must be a FluidBoundaryCondition. 

356 # 

357 # Returns 

358 # ------- 

359 # np.float64 : 

360 # The resistance. 

361 # """ 

362 # if isinstance(node, Node): 

363 # assert isinstance(node.material, Fluid) 

364 # 

365 # else: 

366 # assert isinstance(node, FluidBoundaryCondition) 

367 

368 

369class HeatRadiation(Resistor): 

370 def __init__( 

371 self, 

372 capacitor_1: TemperatureNode, 

373 capacitor_2: TemperatureNode, 

374 view_factor_12, 

375 view_factor_21, 

376 emission_coefficient_1=1, 

377 emission_coefficient_2=1, 

378 area_1=None, 

379 area_2=None, 

380 ): 

381 """ 

382 Calculates the resistance caused by radiative heat transfer between two capacitors. 

383 

384 This Resistor connects itself to the passed Capacitors. You don't need to use connect() 

385 

386 In numeric terms, this is redundant work because we first calculate the heat flux to determine the 

387 resistance, which we then use in the solver to calculate the heat flux again. However, when analyzing the 

388 system from a control engineering perspective, radiative heat transfer should be included in the system 

389 matrix (A-matrix) rather than the input matrix (B-matrix). This is different from what would happen if we 

390 input the heat flux directly to the capacitor as an internal heat source. 

391 

392 The used equation can be found as equation (17b) in 

393 Stephan et al.: 'VDI-Wärmeatlas', 2019, Springer Vieweg. 

394 Chapter 'K1 Wärmestrahlung technischer Oberflächen - 2.2 Sichtfaktoren/Einstrahlzahlen' 

395 DOI: https://doi.org/10.1007/978-3-662-52989-8 

396 

397 Parameters 

398 ---------- 

399 capacitor_1 : TemperatureNode 

400 The first capacitor or BoundaryCondition. 

401 capacitor_2 : TemperatureNode 

402 The second capacitor or BoundaryCondition. 

403 view_factor_12 : float 

404 View factor from area of capacitor 1 to area of capacitor 2. 

405 view_factor_21 : float 

406 View factor from area of capacitor 2 to area of capacitor 1. 

407 emission_coefficient_1 : float, default=1 

408 Emission coefficient from area of capacitor 1. 

409 The same proportion is absorbed (Kirchhoff's law of thermal radiation) 

410 emission_coefficient_2 : float, default=1 

411 Emission coefficient from area of capacitor 2. 

412 The same proportion is absorbed (Kirchhoff's law of thermal radiation) 

413 area_1 : float, optional 

414 The area of the first capacitor that participates in the heat exchange (in m^2). 

415 area_2 : float, optional 

416 The area of the second capacitor that participates in the heat exchange (in m^2). 

417 """ 

418 assert area_1 is not None or area_2 is not None, "At least one area must be provided." 

419 if area_1 is None: 

420 area = area_2 

421 # switch other parameters, because the equation is always assuming area_1 (and other parameters accordingly) 

422 capacitor_1, capacitor_2, view_factor_12, view_factor_21 = ( 

423 capacitor_2, 

424 capacitor_1, 

425 view_factor_21, 

426 view_factor_12, 

427 ) 

428 # this wouldn't be necessary, but it should match the other values 

429 emission_coefficient_1, emission_coefficient_2 = emission_coefficient_2, emission_coefficient_1 

430 else: 

431 area = area_1 

432 

433 self.capacitor_1 = capacitor_1 

434 self.capacitor_2 = capacitor_2 

435 self.temperature_symbol_1 = capacitor_1.temperature_symbol 

436 self.temperature_symbol_2 = capacitor_2.temperature_symbol 

437 self.area = area 

438 self.emission_coefficient_1 = emission_coefficient_1 

439 self.emission_coefficient_2 = emission_coefficient_2 

440 self.view_factor_12 = view_factor_12 

441 self.view_factor_21 = view_factor_21 

442 

443 self.heat_flux: Expr = ( 

444 Stefan_Boltzmann 

445 * self.area 

446 * self.emission_coefficient_1 

447 * self.emission_coefficient_2 

448 * view_factor_12 

449 / ( 

450 1 

451 - (1 - self.emission_coefficient_1) 

452 * (1 - self.emission_coefficient_2) 

453 * view_factor_12 

454 * view_factor_21 

455 ) 

456 * (self.temperature_symbol_1**4 - self.temperature_symbol_2**4) 

457 ) 

458 resistance: Expr = (self.temperature_symbol_1 - self.temperature_symbol_2) / self.heat_flux 

459 super().__init__(resistance=resistance) 

460 

461 # connect self to both capacitors 

462 self.double_connect(self.capacitor_1, self.capacitor_2) 

463 

464 

465def resistance_bc_heat_transfer(bc: BoundaryCondition, node: Node): 

466 """ 

467 Returns the resistance of a heat transfer between FluidBC-Solid Node or SolidBC-Fluid Node. 

468 

469 Parameters 

470 ---------- 

471 bc : BoundaryCondition 

472 node : Node 

473 

474 Returns 

475 ------- 

476 np.float64 

477 """ 

478 effective_area = node.get_area(bc) 

479 assert is_set(bc.heat_transfer_coefficient), "FluidBoundaryCondition has to get a heat transfer coefficient." 

480 return 1 / (bc.heat_transfer_coefficient * effective_area) 

481 

482 

483def resistance_channel_heat_transfer(channel_node: ChannelNode, node: Node): 

484 effective_area = channel_node.get_pipe_area(node) 

485 assert isinstance(channel_node.material, Fluid) 

486 

487 # calculate the alpha using Gnielinski 

488 alpha = alpha_forced_convection_in_pipe(channel_node.diameter, channel_node.velocity, channel_node.material) 

489 return 1 / (alpha * effective_area)