Theory manual

The engineering basis of the Thermal Architect engine: what it models, how it solves those models, the correlations and property data it uses, where they apply, what has been verified and validated, and what it does not attempt to model.

1. Purpose and scope
2. Modelling philosophy
3. Model hierarchy and fidelity
4. Conventions and terminology
5. Thermal network formulation
6. Steady-state solution
7. Transient solution
8. Solver outcomes and what they mean
9. Thermal resistance elements
10. Convection
11. Contact conductance
12. Radiation
13. Heat pipes and vapour chambers
14. Printed circuit boards
15. Fluid networks
16. Thermal-fluid coupling
17. Material and fluid properties
18. Correlation transparency
19. Correlation register
20. Studies: sweeps, optimisation, envelope
21. Verification
22. Validation
23. Uncertainty and sensitivity
24. Calibration and hardware correlation
25. How much should I trust this result?
26. Assumptions
27. Known limitations
28. When to use CFD, FEA or testing instead
29. Modelling guidance

1. Purpose and scope

Thermal Architect is a system-level, reduced-order engineering tool for thermal and thermal-fluid design. It solves:

• lumped-parameter thermal resistance–capacitance (RC) networks, steady and transient;
• printed circuit boards as two-dimensional finite-volume conduction meshes with effective (homogenised) properties, assembled into the same RC network;
• steady, incompressible, one-dimensional fluid networks (pipes, ducts, fittings, pumps, fans, heat exchangers), coupled to the thermal network through component walls.

It answers questions of the form "how hot does this part run, which path limits it, and how sensitive is that to the inputs" for a design described by its conduction paths, convection and radiation to its surroundings, and the coolant that carries heat away. Heat-transfer coefficients come from published correlations for idealised geometries, or are entered by the user. They are never computed from a resolved flow field.

This manual is a theory and engineering-reference document. It is not a user guide. It is written so that an experienced thermal engineer can reproduce or independently check a calculation, and can decide when a result should and should not be relied on.

The manual must not be read as a statement that the engine is validated for any particular application. Section 22 sets out what has and has not been compared with physical measurement.

2. Modelling philosophy

• Reduced order by design. A body is represented by as few isothermal nodes as describe its important gradients. Physics that a lumped network cannot represent, such as flow fields, local hot spots inside a node, or geometric visibility, is either supplied as an input (a coefficient, a view factor) or not modelled. It is not approximated silently.
• Every number has a stated provenance. Each correlation returns its name, source, dimensionless groups and, where the engine checks it, whether those groups lie inside the correlation's range. The checked/unchecked status of each range is listed in the register.
• Refuse rather than guess. Inputs the engine cannot solve meaningfully are refused with an error: compressors and compressible duct flow, unresolved board laminates on import, structurally floating networks, unknown units in an IPC-2581 file.
• Numerical convergence is necessary but not sufficient. A converged answer is a correct solution of the model as built. It says nothing about whether the model represents the hardware (section 25).

3. Model hierarchy and fidelity

The engine's physics sits at four levels of abstraction. Each level inherits the assumptions of the one below it.

Thermal Architect is a system-level, reduced-order thermal and thermal-fluid engineering modelling tool. It is not, and does not attempt to be:

• computational fluid dynamics (CFD) of air or liquid, laminar or turbulent;
• three-dimensional conjugate heat transfer;
• detailed finite-element analysis (thermal or structural; it computes no stress or deformation);
• compressible-flow analysis (choking, shocks, compression heating, pressure waves);
• geometry-based or ray-traced radiation analysis (view factors are entered or taken from closed forms; visibility and obstruction are not resolved);
• electromagnetic or electrical analysis of a PCB (nets and copper current are not read; Joule heating in traces is not computed);
• two-phase transient simulation of heat pipes (start-up, dry-out recovery);
• a qualification or certification tool, or a substitute for hardware test.

These exclusions describe the tool's intended level of abstraction. They are not defects. Section 28 describes when a different tool answers the question better.

4. Conventions and terminology

Units

Inputs and results are SI. Temperatures are entered and reported in degrees Celsius and converted to kelvin wherever a relation needs absolute temperature (radiation, gas density, property evaluation, rarefaction). Fluid-network pressures are absolute, in kilopascals at the interface.

Sign conventions

Heat flow through a link is reported positive from its first (source) end to its second (target) end. Mass flow through a fluid component is positive from its source node to its target node. Conductive links are bidirectional: heat flows from the hotter end, and the stored order fixes only the sign of the reported value. Advection links are one-way (section 9). Two-phase devices take the first end as the evaporator by convention and reverse roles if heat flows the other way.

Terminology

These terms are distinct in this manual and should not be used interchangeably.

Verification and validation

Verification asks whether the equations are solved correctly. It uses analytical solutions, conservation checks, limiting cases, symmetry, mesh and time-step convergence, and independent re-derivation of the implemented equations. Validation asks whether the equations represent physical reality adequately for an intended use. It needs comparison with measurement. This manual uses the two words only in these senses (sections 21 and 22).

Traceability

Every response and report carries the product version and the engine version. The engine version identifies the physics that produced a number: two results computed by different engine versions can differ where the physics between them changed.

5. Thermal network formulation

Nodes

• A component node is an isothermal lumped mass with thermal capacitance C_i [J/K] and a heat load Q_i(t) [W]. Its temperature is an unknown. A node with C = 0 is algebraic: in a transient it follows the network with no lag.
• A boundary node has a prescribed temperature (ambient air, a cold wall, a controlled plate). It supplies or absorbs whatever heat the network requires and is equivalent to infinite capacitance.

Energy balance

Conservation of energy at each component node i:

with j over linked component nodes, b over linked boundary nodes, and s over external sinks: a coolant stream presented as a conductance to its inlet temperature (section 16). Advection links contribute a one-sided term (section 9). In matrix form, with G the assembled conductance matrix and q the loads plus boundary and sink terms:

G is a weighted graph Laplacian of the links plus a diagonal of boundary and sink conductances. For networks without advection it is symmetric. With non-negative conductances and at least one path to a boundary it is a non-singular M-matrix (positive diagonal, non-positive off-diagonal entries, diagonally dominant). This is why a correctly anchored linear network has a unique solution, and why the discrete solution obeys a maximum principle: with no loads, no node is hotter than the hottest boundary. Advection links make G non-symmetric but keep it diagonally dominant.

Capacitance

A node's capacitance is entered directly, or derived as C = m c_p(T) from a library material and a mass. In the second case it follows the material's specific heat with temperature (section 7 gives when it is evaluated). Heat pipes and vapour chambers add their envelope and working-fluid capacitance to their end nodes in transients.

Heat loads

A load may be constant, a step, a ramp, a periodic pulse train, a sinusoid or a table. A steady solve freezes every load at t = 0 by default and reports each load's cycle average alongside. It answers "where would this settle if the load held at its t = 0 value", not "what is the cycle-averaged temperature".

Floating groups and anchoring

Before a linear solve the assembled matrix is checked structurally. A row whose diagonal exceeds the sum of its off-diagonal magnitudes (it has a conductance to a boundary, a sink or, in a transient, its own capacitance) anchors everything connected to it. Unknowns not reached from an anchor form a floating group, whose temperature is undetermined, and the solve is refused with a message. The check is applied to the matrix as assembled at that sweep or step, so it also catches a path that a temperature-dependent conductance has reduced to zero.

Numerical conditioning

Networks that mix very large and very small conductances are ill-conditioned. Examples are a near-perfect "bond" link of 10⁶ W/K beside a 10⁻³ W/K insulation path, or a large-capacitance node with a very short time step. Direct solvers then lose significant digits in proportion to the condition number. Iterative solvers (below) meet their residual tolerance but may carry a larger error in temperature. No condition-number estimate is reported. Where a path is effectively perfect, merge the nodes rather than joining them with an extreme conductance.

6. Steady-state solution

Linear systems

The matrix is always assembled in sparse form. The linear solver depends on size:

• up to 250 unknowns: dense Gaussian elimination with partial pivoting;
• above 250, where SciPy is installed: sparse direct LU factorisation (SuperLU);
• above 250, without SciPy: Jacobi-preconditioned conjugate gradients for symmetric matrices, or Jacobi-preconditioned BiCGSTAB where advection links make the matrix non-symmetric. Iteration stops at a relative residual of 10⁻¹⁰ (measured against the right-hand side), within an iteration cap. A small system that fails to converge falls back to dense elimination; a large one (above 1,500 unknowns) raises an error.

Direct solves are accurate to rounding error amplified by the matrix condition number. The iterative tolerance bounds the residual, not the temperature error. On well-conditioned networks the three paths agree to well below engineering significance.

Non-linear networks: Picard iteration

The conductance matrix depends on temperature when the network contains radiation, contact, correlation-driven convection, heat sinks, boiling, heat pipes, vapour chambers, thin-film stacks, advection with temperature-dependent c_p, or any link or capacitance backed by a temperature-dependent material. Such networks are solved by Picard (successive-substitution) sweeps on the conductances:

1. Evaluate every link conductance G(T^k) at the current estimate. For radiation this is the secant form in section 12; for two-phase devices it is the secant Q/ΔT.
2. Solve the linear system G(T^k) T* = q for the sweep's target T*.
3. Update with relaxation factor ω: T^k+1 = T^k + ω(T* − T^k).
4. Stop when max_i|T*_i − T^k_i| < 10⁻⁶ K (default), or after 60 sweeps (default).

A network with no temperature-dependent element is solved in a single pass.

Convergence safeguards, as implemented

• Global under-relaxation. ω starts at 1. When the step norm grows from one sweep to the next, ω halves, to a floor of 0.1. When sweeps contract again (step below 0.7 of the previous), it recovers by a factor of 1.5, up to 1.
• Conductance damping. For correlation-driven convection and heat-sink links, each sweep's conductance is a weighted geometric mean of the new and previous values. The weight is 1/(1 + m), where m is the local exponent of G ∝ |ΔT|^m estimated from the last two evaluations and clipped to [0, 3]. This removes the two-cycle that plain substitution produces for natural convection (m ≈ 0.25 to 1). The same damping is applied to radiation links once sweeps alternate or diverge. At a fixed point the damped and undamped conductances coincide, so the converged answer is unaffected.
• Aitken Δ² extrapolation. When successive step vectors have a steady ratio r (changing by less than 5 % between sweeps, with 0.05 < |r| < 0.95), the remaining geometric series is summed in one step, T += δ·r/(1 − r). Later sweeps must still meet the tolerance, so an extrapolation that overshoots costs extra sweeps but does not change the converged answer.
• Boiling links use the previous sweep's heat flux to evaluate the boiling curve (a lagged evaluation).

What the convergence test measures

The stopping criterion is the change in temperature between successive sweeps. It is not the energy residual. For a geometrically contracting iteration with ratio r, the remaining error is roughly |r/(1 − r)| times the last step. At r = 0.9 that is nine times the step, so a 10⁻⁶ K step tolerance still implies a small absolute error. The step is measured to the undamped target, so under-relaxation cannot make the iteration look converged when it is not. Separately, every result carries a node energy balance: the net heat into each node from its load, its links and any external sinks. At a converged steady state the component-node entries should be near zero and the boundary entries carry the heat leaving the model. This balance is the conservation check, and a report should be read with it.

Property evaluation in links

A link's material properties are evaluated at the arithmetic mean of its two end temperatures. For conduction with a conductivity k(T) the exact one-dimensional result uses the integral mean (1/ΔT)∫k dT. The two agree to second order in ΔT. The error is negligible when k varies little across the link and grows with the curvature of k(T) and the size of ΔT. If a link carries a large temperature drop through a strongly temperature-dependent material, divide it into several links. Outside a material's stated temperature range the property is held at its value at the range limit, not extrapolated, and the use is flagged against the link (section 17).

Termination by temperature threshold

If any node temperature exceeds 5,000 °C on three consecutive sweeps, the iteration stops. This is a numerical and physical-validity safeguard, not a proof that no steady state exists. Where a two-phase device is at capacity and the heat beyond it has no other path out, the energy balance genuinely has no bounded solution, and the result says "No steady state" and names the device. Otherwise the result says the solve stopped, and lists the possible causes: a load with no adequate path out, a conductivity or correlation driven far outside its range, an unphysical input (a wrong unit or a missing link), or an iteration that diverged although a bounded solution exists.

Non-convergence

A solve that does not meet the tolerance within the sweep limit is returned marked not converged with the last step size and a note. The result is the last iterate. Read it as an indication of where the solution is, not as the solution. Section 8 distinguishes the possible causes.

Uniqueness

A linear network has exactly one solution. A non-linear network may in principle have more than one, for example with boiling curves, heat-pipe saturation, or strongly non-monotonic conductances. Picard iteration finds one solution, which depends on the starting temperatures. The engine does not search for other solutions. When the outcome matters, solve again from different initial temperatures.

7. Transient solution

Scheme

Transients are integrated by backward (implicit) Euler. For step n → n+1 over [t_n, t_n+1], Δt = t_n+1 − t_n:

• Time step. Fixed and user-chosen. The run is ⌈duration/Δt⌉ steps, and the last step is shortened to end exactly at the requested duration. At most 200,000 steps are allowed. There is no adaptive step control and no local error estimate.
• Loads. q̄ is each load's exact time-average over the step, so the energy delivered matches the profile even for a pulse shorter than a step. The shape of such a pulse's response is still smeared over the step.
• Conductances and capacitances of non-linear elements are evaluated once per step at the start-of-step temperatures Tⁿ, and the step is solved once. There is no inner iteration to make them consistent with Tⁿ⁺¹. The scheme is therefore a linearly implicit (lagged-coefficient) backward Euler method. It is fully implicit only for linear networks.
• Linear networks (fixed conductances and capacitances): C/Δt + G is factorised once and each step is a back-substitution.
• Initial state. The run starts either from the nodes' entered initial temperatures or from the steady solution at t = 0.
• Two-phase devices. The envelope and charge capacitance of heat pipes and vapour chambers is computed once, at the initial temperatures, and added to the end nodes. The device's capacity limit is applied every step at the lagged temperatures. Start-up, freezing and dry-out recovery are not modelled (section 13).
• Coolant loops are re-solved quasi-steadily during a transient: at the start, and at the start of any step at which a ported wall has moved more than 0.01 K since the last solve. The coolant itself stores no heat (section 16).

Stability is not accuracy

For a linear network, backward Euler is A-stable (indeed L-stable). No step size makes the computed solution grow without bound, and stiff modes decay rather than oscillate. For the lagged non-linear scheme, each step solves an M-matrix system, so temperatures stay bounded by the step's sources and boundaries, and no step-size instability has been observed. Strictly, though, A-stability is a property of the linear problem.

Stability does not imply accuracy. Backward Euler is first-order: the global error is proportional to Δt. A step much larger than the relevant time constants produces a smooth, stable, plausible-looking curve that is wrong. It under-predicts peaks, smooths fast transients, lags the response, and misrepresents non-linear behaviour that occurs within a step, such as a heat pipe reaching its limit or radiation rising steeply. Validation case A07 measures the order (observed 0.99) on a single-node problem. A06 shows an error of about 0.1 % of the excursion at the chosen step.

Choosing the time step

• Identify the fastest time constant that matters to the quantity of interest, τ ≈ C_i/ΣG_ij for the node concerned, and the fastest change in any load or boundary (pulse width, ramp duration, sinusoid period).
• As a starting point, use Δt ≤ τ/10 for the fastest relevant node and resolve each pulse with at least 8 steps. For first-order backward Euler the peak error after a step change is roughly Δt/(2τ) of the excursion.
• Nodes much faster than the quantity of interest (small capacitance, large conductance) may be left under-resolved: backward Euler damps them correctly to their quasi-steady value.
• For any result that matters, run a time-step convergence check: halve Δt and confirm the quantity of interest changes by less than the accuracy needed. The difference between the two runs approximately estimates the remaining error of the finer run.

The engine warns when a pulse is shorter than four steps, when a sinusoid has fewer than 12 steps per period, and when a step load or pulse period falls outside the run. It also warns when the step exceeds 0.3 of the local time constant C_i/ΣG of any node holding at least 5 % of the network's total capacitance, measured at the first step (an error of roughly 15 % in the lag after a change). With constant loads, a node that settles within the first 2 % of the run is not warned, because only its start-up is affected. The local time constant is an estimate of the right scale for the step, not an eigenvalue analysis, and the warning is not an error estimate.

Reported history

The solver integrates every step but returns at most 400 samples per series, evenly thinned, with the first and last kept. A plotted history can therefore pass under a short peak. The maximum, minimum (each with its time) and time-mean of every node temperature and link heat flow, and of the hottest solved node, are computed from every step and returned separately as the run's statistics. Study metrics that reduce a transient (max, min, range, mean) use these full-resolution figures.

Energy accounting

The energy delivered by each load is accumulated from the step-mean loads actually applied, and the mean power reported is that energy over the duration. Validation case A09 checks first-law closure over a run to 0.2 %.

8. Solver outcomes and what they mean

These conditions are different and must not be confused:

9. Thermal resistance elements

Each link resolves to a conductance G [W/K] at the current temperatures:

Conductance multiplier. Any link may carry a multiplier on its conductance. Calibration and sensitivity studies use it to scale a path without editing the underlying physics. A multiplier fitted by calibration is an empirical correction. It should be reported as such, not as a corrected correlation.

Spreading and constriction. These are not automatic. A small source on a larger plate represented by a single kA/L link omits the spreading resistance and is optimistic. Add a spreading link, use the PCB mesh, or divide the plate into nodes.

10. Convection

Convection is represented in one of three ways:

1. Entered coefficient (G = hA). Its accuracy is the user's.
2. Correlation link. A published correlation for an idealised geometry: natural convection on plates, cylinders and spheres; enclosures and cavities; forced external flow; tube banks; mixed convection; internal flow including microchannels; and impinging jets. Heat sinks have their own models.
3. Fluid-network wall (section 15). Internal-flow Nusselt correlations driven by the solved flow.

For a correlation link, the first end is the surface and the second the fluid (free stream, bulk, or for an enclosure the opposite wall). Fluid properties are taken at the film temperature (T_s + T_∞)/2 for external, natural, mixed and jet cases; at the bulk temperature for internal flow and microchannels; and at the free-stream temperature for the Zukauskas and Whitaker forms, which carry their own wall-viscosity correction.

What a correlation coefficient is. It is a surface average for the idealised geometry, flow and boundary condition of the original experiments. These are typically an isolated, isothermal or uniform-flux surface with a known approach flow. Published scatter for widely used single-phase correlations is commonly ±15 % to ±25 % inside their range, and larger for natural convection in real enclosures, mixed convection, jets and heat sinks with bypass. In a real product the dominant uncertainty is usually the applicability of the idealisation: how much air actually passes the surface, from what direction, and at what temperature. That uncertainty is not reduced by a more accurate correlation. Treat convection coefficients as uncertain inputs and bracket them (section 23).

Radiation from convecting surfaces is separate. A natural-convection surface in air often loses a comparable amount by radiation, and that requires its own link.

11. Contact conductance

A pressed, nominally flat joint conducts in parallel through the contacting asperities and through the gas in the interstitial gaps:

σ is the RMS combined roughness, m the combined asperity slope, k_s the harmonic-mean conductivity of the two solids, H_c the contact microhardness of the softer surface, E' the effective elastic modulus, Y the mean-plane separation and M the gas rarefaction parameter. The gas conductivity comes from the material database, which is CoolProp-backed for air, helium and nitrogen, at the joint temperature. Gas pressure enters only through M.

What the numbers are worth

Contact conductance is very often one of the largest uncertainties in a thermal model, and these correlations are order-of-magnitude estimators, not hardware predictions:

• Published correlations for the same nominal joint commonly differ from each other by a factor of two, and from measurement by more. Scatter in measured data for nominally identical joints is itself large.
• The dominant unknown is the contact microhardness of the prepared surface. For shallow indentation it typically runs two to four times the bulk Vickers hardness that the material library holds, and it depends on surface preparation and work hardening. When no microhardness is entered the bulk value is used. Because h_s ∝ H_c^−0.95, the joint then reads high: an optimistic estimate, and the result says so.
• The models assume nominally flat surfaces with Gaussian, isotropic roughness. They do not represent waviness, flatness error, bolt-pressure non-uniformity, oxide films, coatings or plating, contamination, creep or relaxation under thermal cycling, loading hysteresis, or interstitial greases and pads (use a TIM link for those).
• Nominal pressure P is assumed uniform over the apparent area. Around bolted joints the real pressure is concentrated near the fasteners.

Relative pressure P/H_c above 0.2, combined roughness outside 0.1 to 10 µm and slopes outside 0.02 to 0.4 are warned. Validation case C01 is the only comparison with measured joint conductance (stainless steel 304 in vacuum). It places both the bulk-hardness and the 3×-microhardness results inside a published band that spans a factor of ten. This rules out an order-of-magnitude error and cannot tell the two estimates apart. Treat a computed contact conductance as a range in a sensitivity or Monte Carlo study. Where the joint governs the design, measure it.

12. Radiation

Physics and its representation

Net exchange between grey, diffuse, opaque surfaces at absolute temperatures T₁, T₂ is Q = σ SS (T₁⁴ − T₂⁴), with SS [m²] the total exchange area for the pair. The engine writes this as a conductance using the algebraic factorisation

This secant form is not an approximation at convergence. At the converged temperatures, G_rad(T₁ − T₂) equals the Stefan–Boltzmann exchange exactly, because the factor is re-evaluated each sweep. It is lagged within the iteration and, in a transient, by one time step. The modelling approximations lie in SS, in the grey-diffuse assumption, and in the isothermal-surface assumption.

Radiation links: three forms of SS

• Single emissivity and view factor (default): SS = ε A F. This is the exact result for a small convex grey body in a large enclosure (F = 1, the enclosure effectively black). For two finite grey surfaces it omits inter-reflection and treats the second surface as black, so it over-predicts exchange. For two infinite parallel plates at ε = 0.6 the over-prediction is a factor of 1.4 (validation case L01 documents the effective-emissivity workaround).
• Two-surface grey exchange: SS from the closed-form two-surface enclosure result (Incropera Eq. 13.23) with each surface's own emissivity and area. It is exact for an isolated two-surface grey diffuse enclosure.
• Total exchange area entered directly, as generated by an enclosure (below).

Radiation enclosures

For N opaque, grey, diffuse, isothermal surfaces the radiosity equations

are solved once per surface with unit emissive power on that surface alone. This gives Hottel's total exchange areas SS_ij, so that q_i = Σ_j σ SS_ij (T_i⁴ − T_j⁴). The enclosure is then expanded into pairwise radiation links carrying those exchange areas. The SS matrix depends only on geometry and emissivity, which are constant, so it is computed once. A re-radiating (adiabatic) surface stays in the solve with no net heat. An opening is a black surface at the surroundings' temperature.

View-factor closure. Each row of the entered view-factor matrix must sum to one within 5 %. A row further off is refused unless the enclosure names a surroundings node to take the remainder. Rows within the tolerance are rebalanced by Sinkhorn–Knopp scaling so that reciprocity (A_iF_ij = A_jF_ji) and summation both hold. The rescaling hides view-factor errors of up to 5 %. It does not correct them.

View factors

View factors are entered, or computed from closed-form catalogue results (Incropera Tables 13.1–13.2, Howell's catalogue, Hottel's crossed strings), or computed numerically for two planar rectangles in any relative position. The numerical method evaluates the double-area integral by Gauss–Legendre quadrature with panel doubling and Richardson extrapolation. It is a numerical integration of an exact expression, and it does not consider obstruction by third surfaces. The engine has no geometric model of the product. It does not determine which surfaces see each other, so visibility, shadowing and partial obstruction are the user's responsibility.

Not modelled

Spectral (non-grey) and directional (specular) behaviour; participating media (gases, smoke); solar and environmental (planetary, albedo) loads, except as entered boundary loads; temperature gradients within a radiating surface; and the temperature dependence of emissivity.

Sources: Incropera et al., Fundamentals of Heat and Mass Transfer, 6th ed., section 13.3; Howell, Mengüç & Siegel, Thermal Radiation Heat Transfer; Hottel & Sarofim, Radiative Transfer (1967).

13. Heat pipes and vapour chambers

Each device is described at one of three levels. All three present to the network a secant conductance Q/ΔT between their end nodes, subject to a capacity Q_max. Every reported limit and resistance carries a basis (analytical, correlation, manufacturer data, user-supplied or approximation).

1. Level 1, effective resistance: a user-supplied resistance and rated Q_max. A vapour chamber may instead be given an effective conductivity, in which case spreading is estimated with the Lee, Song, Au and Moran (1995) spreading relation.
2. Level 2, manufacturer data: a performance curve of Q against the end-to-end temperature difference, optionally one curve per orientation, and Q_max tabulated against operating temperature and orientation. The resistance is whatever the curve gives at the operating point. Between zero and the first point the curve is taken through the origin. Beyond the last point it is held (the capacity limit then applies). An orientation outside the data uses the nearest curve and says so; it is not extrapolated. An operating point outside the data is reported as extrapolated.
3. Level 3, physics-based limits: geometry, working fluid and wick. Saturation properties come from CoolProp. Each transport limit is a correlation-based estimate, and the smallest governs:
   • Capillary: available capillary pressure 2σ cosθ/r_c against the liquid (Darcy) and vapour pressure drops and the axial and normal hydrostatic heads. Tilt is positive when the evaporator is above the condenser (gravity-opposed).
   • Boiling: Chi's nucleation criterion with an assumed nucleation radius of 2.54×10⁻⁷ m. This value is not measured for the device, and the result is sensitive to it.
   • Entrainment: a Weber-number estimate. It is an order-of-magnitude estimate only.
   • Sonic: choked vapour at the evaporator exit.
   • Viscous: the low-temperature (start-up) regime.
   • Condenser: the link's own condenser film and sink.
   • Evaporator heat flux: the smallest of the boiling flux, a datasheet flux, and a typical dry-out flux for the wick structure.

Interpretation of the limits

Each calculated limit is a modelled operating limit: a correlation-based estimate of the heat at which that mechanism would constrain the device. It is not a deterministic failure boundary. Wick permeability, effective pore radius, fill charge, manufacturing variation and the boiling nucleation radius are uncertain. Limit estimates from different references can differ by a large factor, and manufacturers derate published Q_max for these reasons. Where available, a manufacturer's tested Q_max (Level 2) should take precedence over a Level 3 estimate.

Behaviour at and beyond a limit

Below its limit the device conducts through its evaporator, vapour and condenser resistances. When the demanded heat exceeds Q_max, the engine caps the transported heat at Q_max, and any further temperature rise is carried by envelope conduction alone. This keeps Q(ΔT) continuous and monotonic, which the solver needs. Physically, exceeding a capillary or boiling limit generally leads to dry-out, in which the transported heat falls rather than holding at Q_max. Results beyond a limit are therefore optimistic and should be read as "the design exceeds the device's modelled capacity", not as a prediction of the post-limit temperature. A steady state that needs more than the device and the rest of the network can carry ends in the temperature-threshold safeguard (section 6). Below the working fluid's usable range the device is treated as frozen: envelope conduction only.

Transients

In a transient the device's envelope and charge store heat (capacitance fixed at the initial temperatures), and the limit is applied every step. Start-up, frozen start-up, transient dry-out and rewetting, and their time constants are not modelled. A transient through a heat pipe describes a pipe that is already operating.

14. Printed circuit boards

What is discretised

A board is a two-dimensional, two-plate finite-volume model with effective properties. Each face is a plate carrying half the thickness and half the heat capacity, conducting in-plane at the effective in-plane conductivity. The two plates are tied cell by cell through the thickness. For a cell of area A on a board of thickness t:

The board therefore resolves in-plane spreading in two dimensions and one through-thickness resistance per cell. It does not resolve individual copper layers, which layer the copper is on, a single buried plane, conduction along a particular trace, the detailed geometry of a via field, or the edges of a pour. In-plane conduction is shared equally between the two faces. This idealisation can misplace heat between the faces of a board whose copper is concentrated on one side.

Effective properties

For a hand-entered stack-up, in-plane conductivity is copper and laminate in parallel (thickness-weighted), and through-plane is the two in series. A stack-up value is therefore an effective property of the whole board. It assumes every copper layer is continuous over the cell, which over-states in-plane spreading for sparse or split layers. Via arrays set the local through-plane conductivity from barrel plating, fill and the remaining laminate in parallel (validation A27; B07 compares a single via with a published worked calculation). Copper pours add sheet conductance to one face (A28). A conformal coat adds its thickness over its conductivity in series with the film, except over keepouts (A29).

Components reach the board through their junction-to-board resistance θ_JB, spread over the cells under the footprint by area, and where given through θ_JC into whatever is mounted on the part. The two paths act in parallel from the junction (A20). θ_JB and θ_JC are JEDEC-defined figures measured in standard environments (JESD51). Using them in a different environment is a standard approximation and is uncertain, typically by tens of percent. They are not compact thermal models.

Mesh

The mesh is rectilinear. The smallest side of the smallest component spans eight cells, and cells grow geometrically away from components by at most a factor of 1.25 per cell, up to a coarse limit set by the board size and its spreading length. The mesh is capped at 4,000 cells. When the cap binds, open areas are coarsened first and component cells last, and the result reports that the mesh was limited. A board's refine factor (1 to 8) divides every cell size, components included, for a mesh-convergence check within the cap. max_cell bounds only the open-area cells.

Mesh convergence evidence, and what it does not show

Mesh convergence and physical validation are different things. The evidence is:

• Analytical verification on idealised boards. A17 (board as a fin against the cosh profile), A18 (point source on a convecting sheet against the K₀ Bessel solution) and A28 (copper pour as a fin) agree with the exact solutions to within about 0.02 K on rises of roughly 1 to 45 K, inside their 2 % to 6 % tolerances.
• Convergence study (engine 4.16.0). Six single-source boards (50 to 150 mm, chips 3 to 10 mm, in-plane k 5 to 40 W/(m·K), h 10 to 50 W/(m²·K), θ_JB 0.5 to 3 K/W) were refined uniformly two- and three-fold. The junction temperature converges at second order (observed order 2.0 to 2.5). Richardson extrapolation put the previous automatic mesh (five cells across a component, growth 1.45) 2 to 5 % of the rise warm. The current setting is 1 to 2 % warm. This study is recorded in the source (pcb.CELLS_ACROSS_COMPONENT). It is not yet a test.
• Grid-refinement case. A19 solves a 100 × 80 mm board with one 10 mm, 1 K/W part on the automatic mesh and again with every cell halved. It checks that the refined mesh is not cell-limited, and requires the Richardson estimate of the automatic mesh's error to be within 2.5 % of the rise.
• A stiff footprint converges slowly. When a part's θ_JB is far below the board's sheet spreading scale 1/(2πk t), its footprint is held at one temperature. The heat then crowds to the footprint's edges and convergence drops to first order. A 10⁻⁴ K/W source on a 20 W/(m·K), 1.6 mm board read 14 % of its rise warm on the automatic mesh. The engine warns when θ_JB is below a tenth of that scale. Re-solve with refine = 2 to check.
• Symmetry (A22, A26) and energy closure (A19, to 10⁻⁶ relative).

The discretisation error of the automatic mesh is in the conservative (warm) direction in every case measured. None of this is evidence that the board model is physically accurate to any stated level. The physical error is dominated by the effective properties, θ_JB, and the film coefficients. These carry much larger uncertainty than the discretisation, and no comparison with a measured board has been made.

Board outlines

A board need not be a rectangle. Its outline is a polygon (arcs taken as chords of at most 10°) with any number of polygonal cut-outs, and the rectangular mesh is laid over its bounding box. Each cell is clipped against the outline and the cut-outs together (Sutherland–Hodgman, even-odd rule). The fraction φ of the cell that is board scales its heat capacity, face convection and through-thickness tie, and each in-plane link is scaled by the fraction of the shared edge that lies on the board. A cell with φ below 10⁻⁴ is dropped. A full-rectangle outline reproduces the rectangular board. The partial-cell treatment is a first-order approximation of the outline: curved edges are stair-stepped at the cell scale.

Imported layouts (KiCad, Altium, IPC-2581)

A KiCad .kicad_pcb, Altium .PcbDoc or IPC-2581 (revisions A to C, .xml/.cvg) design is read into one normalised model: outline and cut-outs, stack-up in physical order, material records, components with their footprints, rotation and side, pads, tracks, copper areas, vias with their spans, and mounting holes. Nothing thermal is taken from the file unless the file states it. No ECAD format carries a thermal conductivity, and IPC-2581 has no property for one in any revision. Each dielectric is therefore mapped by its material name to a laminate datasheet figure or to the material library, or is left unresolved. An unresolved material blocks generation until the user assigns one. Missing copper and dielectric thicknesses, via plating (25 µm) and barrel fill (air) are assumptions, each listed in the import's stack-up table where it can be overridden. Component power, θ_JB, θ_JC and T_j,max are never assumed. A part is included in the model only when the user includes it (or its power reaches the threshold, 0.1 W by default) and it has a power and a θ_JB. Library typical values are applied only on request, at the conservative end of their range.

IPC-2581 specifics. Lengths are converted using the units of the section they appear in (CadHeader, or each dictionary's own). A missing or unknown unit is refused, never guessed. A placement's Xform is applied as scale, counter-clockwise rotation, then mirror in x, then offset, so a mirrored (bottom-side) part's body turns by −rotation in the board frame. Layers are ordered by the stack-up's sequence, not by where elements appear in the file. Negative-polarity copper layers are copper everywhere except their drawn features. Drill-layer holes carry the via spans: plating VIA is a via (through, blind, buried or microvia from its span), PLATED is a component pin's barrel, and NONPLATED is a mechanical hole. Holes of 2 mm or more that are not a pin's are listed as mounting holes.

Component to board. Electrical connectivity is not taken as thermal connectivity, and nets are not read. A part's θ_JB enters the board through the cells its body covers. Where an exposed thermal pad is identified (a pad named EP, PAD or TAB, or one more than four times the area of every other pad of the part), it enters through that pad's area only. This is the concentrated, conservative choice, since θ_JB is measured with the pad soldered. Vias in the pad are counted and reported; their barrels already raise that region's k_z. A BGA's thermal balls are not singled out. Parts can be lumped by region (powers added, θ_JB and θ_JC combined as parallel conductances, the lowest T_j,max kept), and passives can be left out.

Vias. Barrels are never nodes. They can be treated in three ways: spatially, where each region's barrels raise its own k_z as below; as an effective medium, where the same total conductance is spread over the board by area; or ignored. Mounting holes become standoff ports (ordinary board interfaces) only when the user selects one and enters its conductance.

Copper homogenisation. Every copper layer is rasterised (0.1 mm pixels or finer for boards up to about 50 × 50 mm, 4×4 supersampled) into a copper fraction f per pixel. The board is divided into a grid of regions (coarse, medium and fine are 10, 20 and 40 across its longer side). Within a layer, a pixel conducts as copper and resin side by side, k = f k_Cu + (1 − f) k_resin. Across a region, two resistor networks bound the layer's effective conductivity in each direction: rows in series then in parallel, and columns averaged then in series. The region's value is the geometric mean of the two bounds. This choice reproduces the exact result for stripes (where both bounds reduce to the parallel or series limit) and for a two-phase checkerboard (Dykhne's √(k₁k₂)). For general layouts it is an interpolation between rigorous bounds, not an exact value, and its error is bounded by the spread between them. Layers then add in parallel in-plane:

Each plated barrel (vias and plated through-holes) is added in parallel through the thickness, in series with any laminate it does not cross. The strip within 1 mm of the outline is left out of the in-plane bounds, because fabricators pull copper back from the edge, and a series bound across that bare strip would make every edge region an insulator. The regions' k_x, k_y, k_z and volumetric heat capacity become the board's property field. Each mesh cell takes the area-weighted mean of the regions it overlaps, the plates conduct anisotropically (k_x t/2 and k_y t/2), and the mesh is never coarser than a region. A part's θ_JB is spread over the cells its rotated body covers, by area.

15. Fluid networks

The fluid network is the hydraulic analogue of the thermal network. Node pressures play the role of temperatures, component mass flows the role of heat flows, and mass conservation the role of energy conservation. Components are pipes, ducts, fittings, valves, pumps, fans, heat exchangers, chillers and ported walls. Nodes are junctions, and reservoirs at imposed pressure. Each node carries one mixed-mean temperature: the model is one-dimensional, with no profiles across a section.

Flow regime assumptions

• Incompressible. Density is evaluated from the fluid state at each component's inlet and held constant along the component. This is appropriate for liquids, and for gases where velocity is low compared with the speed of sound and the pressure change along the network is a small fraction of the absolute pressure.
• Low-Mach gases. A common engineering rule of thumb is that density changes are below about 5 % for Mach numbers below about 0.3 (roughly 100 m/s in air at room conditions). This rule describes where an incompressible treatment is customarily accepted. It is not a physical limit, and pressure ratio matters as much as velocity: a long, low-velocity duct with a large pressure drop also violates the assumption. The engine does not check Mach number or pressure ratio. The user must check them. Its gas-velocity warning (above 15 m/s) concerns noise in HVAC practice, not compressibility.
• Not represented: compressible flow, choking, shocks, compression and expansion heating of gases (beyond the isenthalpic throttling of liquids, A32), pressure waves, water hammer, surge, and cavitation. Compressors and compressible duct components are refused with an error.
• Steady. The hydraulic solution has no inertia or storage. In a transient it is re-solved quasi-steadily as the walls move, which is valid while the loop's own response is fast against the walls' (section 16).

Hydraulics

Each component's pressure drop, signed with the flow (v|v| so it reverses with the flow), is

• Laminar (Re < 2300): f = C/Re, with C = 64 for a round tube (the Hagen–Poiseuille result, exact for fully developed laminar flow) and the Shah & London value for a rectangular channel. Developing-flow pressure drop in the entrance length is not added.
• Turbulent (Re > 4000): Haaland's (1983) explicit approximation to Colebrook–White, within about 2 % of Colebrook across its range (B05 measures 1.3 % worst case on the tested grid). Colebrook–White itself is a fit to pipe data with an uncertainty commonly quoted around ±15 %.
• Transition (2300 ≤ Re ≤ 4000): f is interpolated linearly between the laminar value at 2300 and the turbulent value at 4000. Transition is intermittent and depends on inlet conditions and disturbances; the interpolation is a modelling convenience, not a physical law. Actual friction in the band can differ substantially, and the engine warns when a component operates there.
• Fittings and valves: tabulated K factors (Idelchik, Crane TP-410; C02 checks the shipped values against published ranges). K is applied as a constant. In reality it depends on Reynolds number at low Re, and interaction between close-coupled fittings is not represented.
• Pumps and fans: from their curves, with the affinity laws for speed and a density correction for fans (A31). Library curves are representative of their class, not a product.
• Hydrostatics: ρ g Δz per component, using that component's density.
• Reverse flow is allowed. Losses reverse sign with the flow. Pump curves are extended for reversed flow and flow past free delivery, and the result flags a machine operating there.

Solution

Junction pressures are solved by Newton–Raphson on nodal mass conservation, with an analytical or finite-difference Jacobian, a step limit, and a backtracking line search (up to 20 halvings). Near-zero flows are linearised about a small fixed reference velocity (10⁻³ m/s), and pump curves about a small fraction of free delivery, so that the Jacobian stays non-singular at zero flow and dead-head. Convergence requires the largest nodal mass imbalance to fall below 10⁻¹⁰ of the largest flow in the network. That is subject to an absolute floor of 10⁻¹² kg/s and a floor set by double-precision resolution of the pressures. At most 100 iterations are made. The final residual and a converged flag are returned.

Energy transport

Once flows are known, enthalpy is carried along them in the flow direction, with perfect mixing at junctions. Heat is added by machine inefficiency, by heat exchangers and chillers, and by ported walls. Properties are re-evaluated at the solved temperatures and pressures, and the hydraulic and energy stages are repeated until flows change by less than 10⁻⁶ of the largest flow, for at most 12 passes (6 for networks with boiling or condensing streams). If they have not settled the result says so.

Wall heat transfer

Heat from a wall at a single temperature T_w into a stream follows the effectiveness–NTU result for a constant-temperature surface:

so a stream cannot leave hotter than its wall, however long the component (A15). The whole wall of a component is assumed to be at one temperature. A long cold plate with a large wall temperature variation should be split into several ported components.

The internal coefficient h comes from:

• Default: fully developed laminar Nu (3.66, constant wall temperature, for a round tube; the aspect-ratio value for rectangular channels) below Re 2300, and Dittus–Boelter above Re 4000, with a linear blend between. Dittus–Boelter is stated for Re ≥ 10⁴ and 0.6 < Pr < 160. Between 4000 and 10⁴ it is extrapolated and warned.
• Opt-in, per component: Hausen's developing-flow laminar value and Gnielinski's correlation above Re 3000, blended between 2300 and 3000. Gnielinski is generally the more accurate of the two turbulent correlations.

Direction of the default's bias. Against Gnielinski (validation B06), Dittus–Boelter reads up to 29 % low for liquids (Pr 3–7) and up to 11 % high for gases (Pr 0.73). The fully developed laminar Nu under-predicts the average coefficient of a developing flow. "Low h" is conservative for a part cooled by the stream (it runs hotter than predicted). It is not conservative where the stream is being cooled, where heat pick-up into the coolant is the quantity of interest, or for gases in the turbulent default. Choose the correlation deliberately.

16. Thermal-fluid coupling

Steady coupled solution

The thermal and fluid networks are solved by partitioned (block Gauss–Seidel, or Picard) iteration:

1. Initial wall temperatures are the fluid temperatures at each ported component's inlet node.
2. The fluid network is solved against the current wall temperatures (including its own property passes).
3. Each ported component is passed to the thermal network as a conductance ṁ c_p ε from its wall node to a sink at the stream's inlet temperature. Components sharing a wall add.
4. The thermal network is solved with those sinks. Its tolerance is loosened while the walls are still moving (between 10⁻³ and 10⁻⁶ K, tracking the coupling movement), and it is warm-started from the previous pass.
5. The wall temperature change is the coupling residual. If its largest entry is below 10⁻⁴ K the coupling has converged, and a final thermal solve is made at full tolerance with the converged sinks.
6. Otherwise the walls are updated with Aitken dynamic relaxation. The factor is recomputed from the last two residuals and bounded to [0.05, 20], so it can under-relax or over-relax. Return to step 2, for at most 40 passes.

A conductance, rather than a heat rate, is passed across the interface because it keeps the alternation stable: the heat into the stream is linear in the wall temperature and is solved implicitly by the thermal side. The editor, the API, the operating envelope and the report all use this routine.

Limits of the coupling. The scheme has no convergence guarantee. Strong feedback between the loops can slow it or make it oscillate. Examples are boiling streams, chillers whose capacity depends on condenser water temperature, and strongly temperature-dependent viscosity. If it does not reach the tolerance in 40 passes, the result is returned with the coupling marked not converged and the final wall movement. Non-linear coupled systems can in principle have more than one consistent solution, and the iteration finds the one reached from its starting point. Energy is conserved at the interface at convergence: the heat leaving each wall through its sink equals the heat the fluid side computes for that component (A23).

Transient analyses with a coolant loop: quasi-steady re-coupling

During a transient the coolant loop is re-solved against the current wall temperatures. This happens at the start, from the initial temperatures (or from the coupled steady state when the run starts from steady), and again at the start of any step at which a ported wall has moved more than 0.01 K since the last solve, up to 300 solves per run. Between solves, each wall sees the loop as a fixed conductance ṁ c_p ε to the coolant inlet temperature. The result reports how many times the loop was re-solved. Consequences:

• A closed loop's coolant warms as the heat it carries comes back round, following the walls with a lag of about the 0.01 K tolerance. On the liquid-cooling example the run's end state is within 0.02 K of the coupled steady state.
• The coolant, its pipes and any reservoir store no heat. Coolant temperatures respond to the walls instantly, so a loop whose own thermal mass matters warms too quickly. Represent such a mass with thermal nodes.
• Flow has no inertia, and changes of flow during the run (a pump stopping, a valve closing) are not represented. The pumps run on their steady curves throughout.
• Because the coupling is lagged and explicit between solves, convergence is not guaranteed for a strongly non-linear loop (boiling streams, chillers). The 0.01 K tolerance keeps the lag small where the walls have thermal mass.

Studies with a coolant loop

Parameter sweeps, optimisation, Monte Carlo, sensitivity and the Engineer Advisor solve the loop once at the base design and hold it fixed while the thermal parameters vary, for steady and transient metrics alike. This is exact for the base design and approximate for changes that move the coolant flow or heat pick-up. The result carries a note. The operating envelope is the exception: it re-solves the full coupled problem at every operating point.

17. Material and fluid properties

Range of validity. Each solid states the temperature range its property models apply over. A material that does not state one defaults to −50 to 300 °C. Outside the range, properties are held at their value at the nearer limit (clamped), because extrapolating a fitted curve can drive a conductivity towards zero and make the solve unstable. The use is flagged against the link.

Property caching. CoolProp calls are cached on state rounded to 0.01 K and 0.001 kPa.

Uncertainty. Alloy temper, ceramic grade and density, polymer filler loading, anisotropy of composites and laminates, and a specific product's pump or fan curve can each move a result by more than the solver's numerical error. Thermal conductivity of commercial alloys of the same nominal designation commonly spans ±10 % to ±20 %. Filled polymers and TIMs span more. Built-in values are suitable for sizing and comparison. Replace them with supplier or measured data before a result is relied on for a final design. Property data are not validated by this engine; CoolProp and the handbooks are the authorities for their own data.

18. Correlation transparency

The register gives each correlation's reference, validity range and whether the engine checks it. This section summarises, by family, the information an engineer needs to judge a result. Uncertainty figures are typical values from the open literature for use inside the stated range. They are not measured accuracies of this implementation, which reproduces the published correlations (section 21).

The uncertainty column is general engineering guidance from the literature, not a statement about this engine.

19. Correlation register

Every correlation the engine offers, read from the engine's own catalogues, with the reference given there and the range over which it applies. The last column says what the engine does about that range:

• warned the engine checks the range. A result used outside it comes back flagged invalid, with a warning naming the quantity. The value is still used in the solve.
• stated the range is from the reference. The engine does not check it, so the user must.
• closed form a closed-form result with no fitted range. The conditions are geometric idealisations (for example "long cylinders", "insulated edges") that the user must confirm hold. "Closed form" means exact for the idealised problem, not for the hardware.

A test in the engine's suite fails if a catalogue gains a correlation with no entry here, so the register is complete by construction.

20. Studies: sweeps, optimisation, envelope

Studies re-solve the same model with changed inputs and add no physics of their own. With a coolant loop, all except the operating envelope hold the loop at its base-design solution (section 16). Transient metrics reduce the full-resolution statistics of each run (section 7).

• Parameter sweep: one or two parameters over a grid (at most 2,000 points), tabulating a scalar metric.
• Optimisation: one parameter at a time. A coarse scan across the stated range (9 points by default) brackets the optimum, then golden-section refinement finds a minimum, maximum or target value within it (at most 200 evaluations), with an optional constraint on a second metric. It finds the best point the scan brackets. If the metric has several optima finer than the scan spacing, it can miss the global one.
• Operating envelope: steady operating points across a range of conditions, each solved fully coupled. Mission profiles (time sequences) can be described but are not solved. There is no transient envelope.
• Engineer Advisor: a heat-path and bottleneck decomposition at the solved, linearised state. Each link's share of a hotspot's rise is R ∂T_i/∂R. For a non-linear network these shares are first-order, at the linearised state. Every recommended change is priced by a full re-solve.
• Monte Carlo and one-at-a-time sensitivity are described in section 23, and calibration in section 24.

21. Verification

Verification shows that the implemented equations are solved correctly. It does not show that they describe the hardware. The evidence is of three kinds.

Unit and regression tests

The repository holds about 1,400 test methods across 54 test files. They cover energy and mass conservation, symmetry, limiting cases, correlation implementations against textbook worked examples, solver paths, import, and regression against stored golden results. Some of them test non-physics parts of the product (accounts, storage, billing). Regression ("golden") tests detect change. They do not show that the stored answer was right.

The benchmark suite, classified

The engine's benchmark suite (shown as "validation cases" in the product) has 46 cases and 148 compared quantities. All pass in the recorded run (engine 4.16.0, 2026-10-01). Under the definitions of section 4, the suite is overwhelmingly verification:

Notes on reading the suite.

• Relative deviations are computed on Celsius values for temperature quantities: (actual − expected)/|expected| with temperatures in °C. Such a deviation depends on the temperature origin, and a percentage tolerance on a Celsius temperature is not a tolerance on the temperature rise. As an example, A18 reports 0.06 % on a temperature of 30.14 °C. That is an error of 0.017 K on a rise of about 10 K above ambient, or about 0.17 % of the rise. Errors should be read in kelvin or as a fraction of the rise.
• A19 uses a board on which a uniform two-fold refinement fits within the 4,000-cell cap, and checks that it does (section 14).
• L01 holds both plates at fixed temperatures and measures the plain link's over-prediction, 1.4, directly. It also checks both exact remedies: the effective emissivity and the two-surface link.
• "Pass" means the stated tolerance or band was met. A wide band (C01 spans a factor of ten) passes a wide range of answers.

The full case list, with every reference value, engine value and deviation, is on the validation page. The product labels these cases "validation". In this manual's terms they are a verification and benchmark suite, plus one coarse experimental bound check.

46 of 46 validation cases pass (engine 4.16.0, run 2026-10-01). 148 quantities are compared. Every case, with its reference and error.

22. Validation

Validation is comparison with physical measurement for an intended use. The current validation evidence is:

• One coarse comparison with measured data (C01, contact conductance of ground stainless steel in vacuum, against a published band spanning a factor of ten).
• Indirect validation through published correlations. Each correlation was fitted to, and validated against, experimental data by its original authors within its stated range. The engine reproduces those correlations (verification), so their published accuracy carries over to their intended use. It does not carry over to a product whose geometry and flow differ from the idealisation.

Not validated:

• any complete thermal model against measured hardware temperatures;
• the PCB model (effective properties, two-plate idealisation, θ_JB coupling) against a measured board;
• the heat-pipe and vapour-chamber models against measured devices;
• the fluid network against measured pressure drop, flow distribution or heat pick-up;
• thermal-fluid coupling against a measured loop;
• transient response against measured transients;
• property data. CoolProp and the handbooks are their own authorities; the engine uses them.

The largest source of error in a real model is almost always whether the model represents the hardware: boundary conditions, convection coefficients, contact resistances, airflow and power. No solver verification addresses this. Calibrating against test data (section 24) and checking the calibrated model against held-out data is how a specific model is validated for a specific use.

23. Uncertainty and sensitivity

Kinds of uncertainty

Monte Carlo

Inputs are sampled from normal (or tolerance), uniform, triangular or lognormal distributions, by simple random or Latin hypercube sampling with a reproducible seed. The default is 500 samples. At most 20,000 are allowed for steady metrics and 400 for transient metrics. Samples are clipped to stated bounds and counted. Inputs are sampled independently; correlation between inputs is not represented. Failed samples are excluded and counted. The result gives the metric's distribution, percentiles, running mean and 95th percentile (to judge whether the sample count was enough), Pearson and Spearman correlations of each input with the output, and with a specification limit the fraction outside it with a 95 % Wilson interval and a normal-fit C_pk.

What Monte Carlo does not do. It propagates the uncertainty of the inputs that are made uncertain, through the model as built. It does not quantify model-form uncertainty, correlation uncertainty that has not been entered as an input, or error from a model that omits a heat path. A narrow Monte Carlo spread is a statement about the inputs, not about the model's fidelity. A normal-fit C_pk assumes the output is normal, which it often is not (check the histogram).

Sensitivity

One-at-a-time sensitivity (a tornado) moves each input down and up by a stated fraction of itself (10 % by default) and re-solves. It is local and ignores interactions. Monte Carlo rank correlation is the global complement.

24. Calibration and hardware correlation

The engine includes a calibration module for the thermal network. It fits chosen parameters to measured data and then assesses what the fit means.

Data

An experiment is a steady or transient run. Its channels map to network nodes: measured temperatures (outputs the model must reproduce), measured or constant power (inputs written into node loads), and measured or constant boundary temperatures. Each temperature channel carries a measurement uncertainty σ; the default is 0.5 K. Preprocessing works on a copy of the data and records every step in an audit log. Each experiment is tagged as a calibration or a validation dataset.

Fit

The residual for each measured point is (model − measured)/σ, weighted so that each channel has equal influence however many samples it logged. Where a parameter has a nominal value and an uncertainty, a prior term is added. Objectives are χ², weighted least squares, RMSE, MAE or maximum error. Optimisers are Levenberg–Marquardt or Nelder–Mead, with multistart (3 starts by default) within parameter bounds. The plain χ² is always reported.

Assessment

• Identifiability: from the Jacobian at the calibrated point, the sensitivity of the readings to each parameter, the parameter covariance (scaled by reduced χ² and by residual autocorrelation, because neighbouring samples of a transient are not independent), the correlations between parameters, and the eigenvectors of the information matrix. Each parameter is graded good, moderate, poor or none by explicit thresholds. Example: two resistances in series observed only at their ends cannot be separated, and the analysis says so.
• Residual analysis: patterns in the residuals that indicate a structural model error rather than noise.
• Validation against independent data: predictions for datasets the fit never saw, and leave-one-experiment-out cross-validation.

What calibration can and cannot do

• Calibration reduces the discrepancy between a model and the measured hardware under the tested conditions. It does not show that the model is physically correct, or that it predicts well outside those conditions.
• A calibrated multiplier absorbs every error in its path: correlation error, contact error, geometric simplification, and errors elsewhere that happen to correlate with it. Report it as an empirical correction factor.
• Fitting more parameters than the data can identify will over-fit: the residuals fall and the predictions get worse. Use the identifiability grades, fix unidentifiable parameters at nominal values, and prefer fewer, physically meaningful parameters.
• To demonstrate predictive capability, validate on data not used in the fit: a different power level, ambient, flow rate or duty cycle. Agreement on the calibration set alone is not evidence of predictive capability.
• Calibration cannot compensate for an unmeasured input that changes between tests, or for thermocouple placement errors. These must be controlled in the test.

25. How much should I trust this result?

A result should be judged on at least five separate dimensions. A model can be excellent on one and poor on another, and a numerically converged result can still be physically unreliable.

A converged result can be physically unreliable when:

• a correlation is used outside its validated range or for a geometry it does not describe;
• contact conductance governs and is uncharacterised;
• material properties are nominal, and the actual alloy, grade or filler is unknown;
• geometry has been simplified past the point where spreading or local gradients are represented;
• airflow distribution (bypass, recirculation, maldistribution) is unknown;
• a body with significant internal gradients is represented by one isothermal node;
• a coolant loop's own thermal mass matters in a transient, or a study varies parameters that move the frozen coolant;
• a two-phase device operates near or beyond its modelled limit;
• there has been no experimental calibration or validation of a similar configuration.

As a rule, the engine's numerical error, once converged, is much smaller than the uncertainty in convection, contact and property inputs. Effort is better spent narrowing those inputs than tightening tolerances.

Recommended result-confidence framework

The following is a recommended interpretation framework. The engine does not currently assign these categories as a single status. It does return the ingredients: converged flags, correlation validity flags, property-range flags, coupling convergence, transition and heat-pipe capacity warnings, and representative-data labels. A future interface could combine them into this status.

A result may carry more than one status, for example VALID numerically and an ESTIMATE in its data. INVALID takes precedence.

26. Assumptions

• Each node is isothermal. Gradients within a body are represented only by dividing it into more nodes, or by a board's mesh.
• Link conductances are evaluated at the link's mean temperature. Film properties are taken at the film, bulk or free-stream temperature as each correlation requires. Properties outside a material's range are held at the range limit.
• Convection coefficients are surface averages from correlations for idealised geometries, or values supplied by the user.
• Radiating surfaces are grey, diffuse, opaque and isothermal. There is no participating medium, view factors take no account of obstruction, and only entered links and enclosures exchange radiation.
• Fluid flow is incompressible, steady and one-dimensional, with one mixed-mean temperature per node and perfect mixing at junctions. Properties are held along a component at its inlet state.
• A wall exchanging heat with a stream is at one temperature along the component.
• In transients the coolant loop is quasi-steady: re-solved as the walls move, with no heat capacity of its own. In studies other than the operating envelope it is held at its base-design solution.
• Transient non-linear conductances and capacitances are lagged by one step.
• Board properties are effective values for the whole stack-up, or per region for an imported layout. In-plane conduction is shared equally by the two faces.
• Contact conductance assumes nominally flat surfaces with Gaussian roughness and uniform pressure, and uses bulk hardness unless a contact microhardness is supplied.
• Heat pipes beyond their limit carry Q_max plus envelope conduction (optimistic relative to dry-out).
• Built-in material, TIM, pump and fan data are nominal or representative.

27. Known limitations

What the engine does not model, what that does to an answer, and what to do about it. The first four are printed on every report.

28. When to use CFD, FEA or testing instead

Reduced-order network models and higher-fidelity methods answer different questions at different levels of abstraction. Neither is inherently better.

Where this engine's approach is appropriate

• early architecture studies and concept selection;
• thermal sizing of heat sinks, cold plates, heat pipes, pumps and fans;
• trade studies and design-space exploration, where many variants must be solved quickly;
• system-level thermal budgets: where the temperature is spent along a path;
• transient response of lumped systems: duty cycles, pulses, soak and warm-up;
• interactions between components, boards and coolant loops;
• sensitivity and uncertainty studies over many inputs;
• coupled thermal and hydraulic system sizing at steady state;
• correlating a system model to test data for prediction at nearby conditions.

Where a higher-fidelity method or test may be required

A common workflow combines them. Size the system here, use CFD or test to determine the convection coefficients and flow split that the network takes as inputs, bring those back, and correlate the network to the prototype test.

29. Modelling guidance

This section gives concise engineering guidance only. Workflows belong in a separate Engineering Guide.

• Node granularity. One node per body is adequate when the body's internal Biot number (hL/k with L its conduction length) is well below 0.1, or when its internal resistance is small compared with the path it sits in. Otherwise subdivide it, and confirm that the quantity of interest stops changing as you do.
• Plain R against a detailed element. Use a plain resistance when the value is known from data, such as a datasheet θ, a measured TIM or a supplier curve. Use a correlation, contact, spreading or heat-sink element when the value must follow geometry, temperature or flow. Do not double-count: a datasheet heat-sink resistance already includes its convection.
• Time step. See section 7. Start at τ/10 of the fastest relevant node, resolve pulses with at least eight steps, then halve the step and compare.
• Convergence. Check the converged flags (thermal, fluid, coupling) and the node energy balance on every result used for a decision. Investigate any not-converged result before using it.
• Correlation validity. Read every warning. For "stated" ranges, check the dimensionless groups by hand. Keep internal flows out of the transition band where possible.
• Default properties. Before relying on a result for a final design, replace library materials, TIMs, pumps and fans with supplier data, and enter or bracket contact microhardness.
• Warnings. A warning means the number was computed but the method was stretched. Quantify the effect with a sensitivity study rather than ignoring it.
• Sensitivity. Run a tornado or Monte Carlo on the convection coefficients, contact and TIM resistances, power and ambient temperature. If the design margin is smaller than the spread, the design is not yet decided by the model.
• Test data. Calibrate on one set of conditions and check the prediction on another. Use calibrated multipliers only near the conditions they were fitted at.