Heat Transfer for CAE Engineers
Conduction, Convection, Radiation & FEM Analysis
Heat management is a critical constraint in virtually every engineering product — from semiconductor chips dissipating hundreds of watts per square centimeter, to jet turbine blades operating above material melting points, to spacecraft facing cryogenic space and 1600°C re-entry simultaneously. Understanding the three heat transfer modes and how they're modeled in FEM/CFD is essential for any CAE engineer.
1. Three Modes of Heat Transfer
I know conduction is heat flowing through a solid, convection involves fluid, and radiation doesn't need a medium. But in a real simulation — say, a heat sink on a power electronics board — which modes are actually important?
For a power electronics heat sink, all three can matter, but in different proportions. Conduction dominates inside the solid metal — from the chip junction through the substrate and into the fins. Convection is the main mechanism from the fin surfaces to the air, and it can be either natural (buoyancy-driven) or forced. Radiation is typically 5–15% of total heat dissipation at moderate temperatures — often neglected in a first analysis, but important to include for accurate results. In contrast, for a spacecraft radiator panel in orbit, radiation is the only mechanism — there's no air to convect with.
| Mode | Mechanism | Governing Law | Key Parameter |
|---|---|---|---|
| Conduction | Molecular vibration / electron diffusion in solids | Fourier's law: $\mathbf{q}'' = -k\nabla T$ | Thermal conductivity $k$ [W/(m·K)] |
| Convection | Bulk fluid motion carrying thermal energy | Newton's cooling: $q'' = h(T_s - T_\infty)$ | Heat transfer coefficient $h$ [W/(m²·K)] |
| Radiation | Electromagnetic wave emission/absorption | Stefan-Boltzmann: $q'' = \varepsilon\sigma(T_s^4 - T_\text{surr}^4)$ | Emissivity ε, view factor F |
2. Heat Conduction Equation
The heat conduction equation — is it similar to the diffusion equation in fluid dynamics? I feel like I've seen this same math in multiple places.
Exactly right — it's the same parabolic PDE structure. The heat equation is the prototypical diffusion equation. The same mathematics governs mass diffusion (Fick's law), neutron diffusion in reactor cores, and moisture diffusion in building materials. Once you understand the heat equation numerically, you understand all of them. The key parameter is the thermal diffusivity α, which controls how fast temperature disturbances propagate through the material.
2.1 General Heat Equation
Expanded in Cartesian coordinates:
For isotropic, constant-property material with no heat generation:
Thermal diffusivity $\alpha$ is the material property that matters for transient problems — it tells you how fast a temperature wave penetrates the material. For copper: $\alpha \approx 1.17 \times 10^{-4}$ m²/s; for stainless steel: $\alpha \approx 3.5 \times 10^{-6}$ m²/s — copper equilibrates ~33× faster.
2.2 Steady-State Conduction (Laplace Equation)
When there's no time dependence and no internal heat generation:
For 1D with heat generation, the solution gives a parabolic temperature profile:
This is the temperature profile in a nuclear fuel rod, in a current-carrying wire, or in the wall of an induction-heated pipe.
2.3 Thermal Resistance Network
For 1D problems, thermal systems can be analyzed exactly like electrical circuits:
For resistances in series (e.g., chip → TIM → heat spreader → air):
In our electronics cooling analysis, the thermal interface material (TIM) between the chip and heat sink seems to have a big impact. How do I model that correctly?
TIM is critical — a poorly applied TIM can add 5–15°C to your junction temperature compared to the theoretical minimum. In FEM, model it as a thin solid layer with the actual k value of the material (thermal grease: 3–8 W/(m·K), phase-change materials: 3–6 W/(m·K), indium foil: 80 W/(m·K)). Alternatively, if the layer is very thin, model it as a contact resistance using a thermal conductance per unit area (W/(m²·K)). The contact resistance R_c'' also includes the effect of microscale surface roughness — even with perfect TIM, there's a residual contact resistance.
2.4 Thermal Conductivity of Engineering Materials
| Material | k [W/(m·K)] | ρ [kg/m³] | c_p [J/(kg·K)] | α [m²/s] |
|---|---|---|---|---|
| Diamond | 2000–2500 | 3515 | 502 | ~1.13×10⁻³ |
| Copper (pure) | 385 | 8960 | 385 | 1.17×10⁻⁴ |
| Aluminum 6061 | 167 | 2700 | 896 | 6.9×10⁻⁵ |
| Steel (carbon) | 50 | 7850 | 490 | 1.3×10⁻⁵ |
| Stainless steel 316 | 16.3 | 8000 | 502 | 4.1×10⁻⁶ |
| Water (20°C) | 0.598 | 998 | 4182 | 1.43×10⁻⁷ |
| Air (20°C) | 0.0257 | 1.204 | 1005 | 2.12×10⁻⁵ |
| Insulation foam | 0.03–0.05 | 30–50 | 1000 | ~1×10⁻⁶ |
| Silicon (chip material) | 150 | 2330 | 700 | 9.2×10⁻⁵ |
3. Convective Heat Transfer
In the software I just enter a heat transfer coefficient h. But where does that number actually come from? My colleague said "just use 10 W/(m²·K) for natural convection in air" — is that reliable?
That "10 W/(m²·K)" rule of thumb is a rough estimate that might be off by a factor of 3 depending on surface orientation, temperature difference, and geometry. For quick screening it's fine. But if your design hinges on convective cooling, use the proper Nusselt number correlations for your geometry — vertical plate, horizontal cylinder, flat plate with forced flow — they give you h as a function of Re, Pr, and geometry. Or better yet, do a conjugate CFD analysis where the solver computes h automatically from the flow field.
3.1 Dimensionless Numbers for Convection
| Dimensionless group | Physical meaning | Typical range |
|---|---|---|
| Nusselt (Nu) | Convective vs. conductive heat transfer; dimensionless h | 1 (pure conduction) to 10⁴+ |
| Prandtl (Pr) | Momentum diffusivity vs. thermal diffusivity | Liquid metals ~0.001; water ~7; oils ~100–1000 |
| Rayleigh (Ra) | Buoyancy vs. diffusion — drives natural convection | Ra < 10⁹: laminar; Ra > 10⁹: turbulent NC |
3.2 Key Nusselt Correlations
Fully developed turbulent pipe flow (Dittus-Boelter, Re > 10,000):
Isothermal flat plate, forced convection (Churchill-Bernstein for all Re):
Vertical plate, natural convection (Churchill-Chu):
Horizontal cylinder, natural convection (Churchill-Bernstein):
3.3 Typical Heat Transfer Coefficient Values
| Configuration | h [W/(m²·K)] |
|---|---|
| Natural convection in air | 2–25 |
| Forced convection in air (fan-cooled) | 25–250 |
| Natural convection in water | 200–1,000 |
| Forced convection in water | 1,000–15,000 |
| Boiling water (pool boiling) | 3,000–60,000 |
| Condensing steam | 5,000–100,000 |
| Liquid metal (NaK) | 10,000–100,000 |
4. Radiation Heat Transfer
I always struggle with radiation. When does it become important enough to model, and how does emissivity actually work in practice?
Good rule of thumb: radiation scales as T⁴, so it gets increasingly dominant at high temperatures. Below 300°C in air, radiation is usually <20% of total heat loss — often safely neglected in a conservative (unsafe) design, but not always. Above 500°C, radiation often dominates over convection. In vacuum (space, furnaces), radiation is essentially the only mechanism. For emissivity: a polished aluminum surface has ε ≈ 0.05 — terrible radiator. The same aluminum with a black anodized coating has ε ≈ 0.85. That's a 17× difference in radiated power. Never use the default ε = 1 without checking your actual surface finish.
4.1 Blackbody and Real Surface Radiation
A perfect blackbody emits the maximum possible radiation at temperature T:
Real surfaces emit less, characterized by emissivity ε (0 ≤ ε ≤ 1):
Net radiation heat transfer between a convex surface 1 and surrounding enclosure at $T_\text{surr}$:
4.2 View Factors
For radiative exchange between two surfaces, the view factor $F_{12}$ is the fraction of radiation leaving surface 1 that directly strikes surface 2:
Key view factor rules:
- Summation: $\sum_j F_{ij} = 1$ — all radiation must go somewhere
- Reciprocity: $A_i F_{ij} = A_j F_{ji}$
- Self-view: $F_{11} = 0$ for flat or convex surfaces
| Configuration | View Factor |
|---|---|
| Infinite parallel plates | $F_{12} = 1$ |
| Concentric cylinders/spheres (inner → outer) | $F_{12} = 1$ |
| Two perpendicular plates sharing an edge, equal width w | $F_{12} = 1 - \frac{1}{\sqrt{2}} \approx 0.293$ |
| Disk to coaxial parallel disk (equal radius, spacing h) | $F_{12} = \frac{1}{2}\!\left[S - \sqrt{S^2 - 4(R_2/R_1)^2}\right]$, $S = 1 + (1+R_2^2)/R_1^2$ |
4.3 Radiation Network Method
For gray diffuse surfaces, radiation can be analyzed using an electrical resistance network:
where $E_{b,i} = \sigma T_i^4$ is the blackbody emissive power and $J_i$ is the radiosity (total radiation leaving surface i, both emitted and reflected).
4.4 Emissivity of Engineering Surfaces
| Surface | Temperature | Emissivity ε |
|---|---|---|
| Polished copper | 20°C | 0.02–0.05 |
| Polished aluminum | 100°C | 0.04–0.06 |
| Black anodized aluminum | 100°C | 0.80–0.95 |
| Oxidized steel | 200°C | 0.70–0.80 |
| Stainless steel, polished | 200°C | 0.10–0.20 |
| Concrete / brick | 20°C | 0.85–0.95 |
| Black paint | 50°C | 0.92–0.97 |
| Bare silicon chip | 100°C | ~0.70 |
5. FEM Thermal Analysis
When I run a steady-state thermal analysis in Abaqus or Ansys Mechanical, what equation is the solver actually solving? Is it similar to structural FEM?
Very similar structure, but simpler — temperature is a scalar (1 DOF per node vs. 3 for structural), so the system is smaller. The conductivity matrix plays the role of the stiffness matrix, and heat flux boundary conditions play the role of force boundary conditions. The big difference is that thermal analysis is often much easier to converge — no contact nonlinearity, no geometric nonlinearity from large deformations. The main source of nonlinearity in thermal problems is temperature-dependent material properties and radiation (T⁴ term).
5.1 Weak Form and FEM Thermal System
Starting from the strong form of the heat equation, the weak form (after integration by parts) is:
Discretizing with shape functions $N_a$, the FEM system is:
where:
- $C_{ab} = \int_\Omega N_a\,\rho c_p\,N_b\,d\Omega$ — thermal capacity (heat capacity) matrix
- $(K_T)_{ab} = \int_\Omega \nabla N_a \cdot k\nabla N_b\,d\Omega$ — thermal conductivity matrix
- $Q_a = \int_\Omega N_a\,\dot{q}\,d\Omega + \int_{\Gamma_q} N_a\,\bar{q}''\,d\Gamma + \int_{\Gamma_c} N_a\,h(T_\infty - T)\,d\Gamma$ — thermal load vector (includes convection boundary condition)
5.2 Convection Boundary Condition in FEM
The convective boundary condition $q'' = h(T - T_\infty)$ contributes both to the stiffness (conductivity) matrix and the load vector:
The assembled system becomes: $(\mathbf{K}_T + \mathbf{H}_c)\mathbf{T} = \mathbf{Q} + \mathbf{f}_c$
5.3 Transient Thermal: Time Integration
The generalized trapezoidal method (θ-method) for $\mathbf{C}\dot{\mathbf{T}} + \mathbf{K}\mathbf{T} = \mathbf{Q}$:
| θ value | Method | Accuracy | Stability |
|---|---|---|---|
| 0 | Forward Euler (explicit) | 1st order | Conditionally stable: Δt < Δt_crit |
| 1/2 | Crank-Nicolson | 2nd order | Unconditionally stable (can oscillate) |
| 2/3 | Galerkin | 2nd order | Unconditionally stable, less oscillation |
| 1 | Backward Euler (implicit) | 1st order | Unconditionally stable, most diffusive |
Practical note: In Abaqus/Ansys, the default transient thermal integration is usually the modified Crank-Nicolson or backward Euler. For problems with sudden thermal shocks (welding, quenching), backward Euler (θ=1) is more robust as it doesn't exhibit spurious oscillations near sharp temperature fronts.
6. Practical Thermal Analysis Tips
I'm getting very high peak temperatures at one corner of my model, but the rest seems reasonable. Could that be a real hot spot or a numerical artifact?
Check these in order: (1) Is there a point heat source applied at that node? Point loads in thermal analysis cause T → ∞ singularities just like stress singularities — distribute the heat over an area instead. (2) Is the mesh very coarse there? Refine and see if the temperature changes. (3) Are there conflicting boundary conditions — a fixed temperature and a heat flux applied to the same surface? That's a setup error. If none of these, it might be a genuine hot spot from internal geometry — look at the local thermal resistance path and ask whether the material and geometry make sense.
6.1 Conjugate Heat Transfer (CHT)
When the fluid flow and solid conduction are both important, conjugate heat transfer analysis couples the fluid energy equation with the solid heat equation at the interface:
CHT avoids the need to specify a heat transfer coefficient h — the solver computes it from the actual flow field. This is more accurate but computationally more expensive. Use CHT when:
- h varies significantly over the surface (complex flow, separation regions)
- The solid temperature distribution affects the flow (buoyancy, transpiration cooling)
- You need to optimize the cooling channel geometry
6.2 Thermal Contact Resistance
At any interface between two solids, surface roughness creates a thermal contact resistance even with mechanical pressure. The interface thermal resistance:
Typical values: metal-metal interface (machined, no TIM): 10⁻⁴–10⁻³ m²·K/W. With thermal grease: 10⁻⁵–10⁻⁴ m²·K/W. In FEM, model as a "thermal contact conductance" (inverse of resistance) applied to the shared surface pair.
6.3 Phase Change: Latent Heat
Melting and solidification involve latent heat that must be accounted for. The enthalpy method avoids the moving boundary problem:
where $L$ is the latent heat [J/kg] and $f_L$ is the liquid fraction (0 in solid, 1 in liquid, 0–1 in mushy zone). This is the basis for solidification modeling in casting simulations.
6.4 Thermal Analysis Checklist
- Material properties: Temperature-dependent k, c_p are important for large ΔT; radiation must use T-dependent emissivity for high-temperature work
- Boundary conditions: Confirm heat loads (W or W/m²?), ambient temperature, and h values are physically realistic
- Steady vs. transient: If the thermal time constant (ρcpL²/k) is much smaller than your load variation period, steady-state is sufficient
- Mesh: Thin walls (PCBs, solder joints) need through-thickness elements; a single layer is often inadequate
- Radiation: Check surface emissivity settings; defaults are often ε = 0 or ε = 1, neither of which may be correct
- Energy balance: Sum of all heat inputs should equal sum of all heat outputs in steady state
- Three modes: conduction ($\mathbf{q}'' = -k\nabla T$), convection ($q'' = h\Delta T$), radiation ($q'' = \varepsilon\sigma T^4$)
- Nusselt correlations (Dittus-Boelter, Churchill-Chu) give h without full CFD
- Radiation scales as T⁴ — critical above 300°C and in vacuum
- FEM thermal system: $\mathbf{C}\dot{\mathbf{T}} + \mathbf{K}\mathbf{T} = \mathbf{Q}$ — analogous to structural dynamics
- Always check energy balance in steady-state: heat in = heat out
- Conjugate heat transfer eliminates the need to specify h — the solver computes it
Related Interactive Tools
- 1D Heat Diffusion Simulator — Animate transient temperature propagation with FTCS explicit scheme
- Fin Efficiency Calculator — Rectangular fin temperature distribution and efficiency η
- Multi-Layer Wall Thermal Resistance — Temperature distribution through composite walls
- Lumped Capacitance Cooling Curves — Biot number check and transient temperature
- Natural Convection Calculator — Churchill-Chu correlation for vertical plate