Heat Transfer · Chapter 5
0%
Heat Transfer · Chapter Five

Numerical Methods in Heat Conduction

Real problems rarely have simple geometries or boundary conditions. Numerical methods replace the governing differential equation with a set of algebraic equations solved at discrete points.

This chapter develops the finite-difference and energy-balance approaches, formulates interior and boundary nodes for one- and two-dimensional steady conduction, and extends the method to transient problems — introducing the explicit and implicit schemes and the stability criterion that constrains the time step.

Learning objectives

  • Translate conduction physics into conservative nodal energy balances.
  • Select mesh spacing, boundary-node control volumes, and a suitable steady or transient scheme.
  • Assemble and solve one- and two-dimensional nodal equations.
  • Apply the explicit stability limit and distinguish stability from accuracy.
  • Verify numerical results using conservation, bounds, symmetry, and mesh refinement.
5–1

Why Numerical Methods

Analytical solutions exist only for simple shapes, constant properties, and idealized boundaries. Most practical cases — irregular geometry, temperature-dependent properties, mixed boundary conditions — demand a numerical approach.

Analytical vs. Numerical

Analytical methods solve the differential equation for a continuous temperature function valid everywhere. Numerical methods instead compute temperatures at a finite set of nodes, trading exactness for the ability to handle nearly any geometry or condition.

Discretization

The medium is divided into a grid of nodal points spaced \(\Delta x\) (and \(\Delta y\)) apart. Each node represents the small volume around it, and the conservation law is written for that element, producing one equation per node.

Physical meaning

A node represents an average over a control volume, not an isolated point. Refining the mesh improves geometric and gradient resolution but increases equation count. Numerical output is credible only after the discrete model, property data, and boundary conditions are audited.

Problem-solving strategy

  1. Recognize steady or transient behavior, dimensionality, generation, and each boundary type.
  2. Sketch control volumes and assign material properties at stated representative temperatures.
  3. Write conservation equations before simplifying coefficients; half and quarter cells matter at boundaries and corners.
  4. Solve the smallest closed system, monitor residuals, and preserve units.
  5. Check global energy balance, maximum-principle bounds, stability, and mesh/time-step independence.
5–2

Finite Difference Formulation

Derivatives are approximated by differences of nodal temperatures using a Taylor-series expansion. The second derivative central to conduction takes a particularly simple three-point form.

Central second-difference approximation
\[ \frac{d^2T}{dx^2}\bigg|_{m} \approx \frac{T_{m-1} - 2T_m + T_{m+1}}{\Delta x^2} \]

Energy-Balance Method

Rather than approximating derivatives, write a steady-state energy balance on the control volume around each node: heat conducted in from all neighbors, plus generation, sums to zero. It gives the same equations as finite differences but with clearer physical meaning and easy handling of boundaries.

Central differences are second-order accurate on a uniform mesh when the temperature field is smooth. Interfaces require conductances based on the series resistance between nodes (often a harmonic-mean conductivity), not a casual arithmetic average. Nonlinear \(k(T)\), radiation, or convection correlations require iteration.

Interior node, 1-D, with generation
\[ T_{m-1} - 2T_m + T_{m+1} + \frac{\dot{g}_m\,\Delta x^2}{k} = 0 \]
5–3

One-Dimensional Steady Heat Conduction

Writing the nodal equation at every interior node, plus special equations at the boundaries, yields a linear system solved simultaneously for all nodal temperatures.

Convection boundary node (energy balance)
\[ h\,A\,(T_\infty - T_m) + k A\,\frac{T_{m-1} - T_m}{\Delta x} + \dot{g}_m\,A\,\frac{\Delta x}{2} = 0 \]
The boundary element is only half the width of an interior element, hence the \(\Delta x/2\) generation term.
EXAMPLE 1Two-Node Temperature Field in a Cooling Rib
Problem. A thin conducting rib is modeled by two interior nodes between boundaries held at \(120^\circ\text{C}\) and \(30^\circ\text{C}\). Uniform spacing, constant \(k\), no generation, and negligible lateral loss apply. Find both nodal temperatures and the normalized discrete residuals.

Strategy

This is steady one-dimensional conduction on a uniform mesh, so each interior balance sets the node equal to the average of its two neighbors.

Givens

  • \(T_0=120^\circ\text{C}\), \(T_3=30^\circ\text{C}\); unknown \(T_1,T_2\).
  • No generation; equal conductances between adjacent nodes.

Physical/model interpretation

Positive heat flows left to right. Constant area and conductivity make the conductance cancel from each energy balance.

Step-by-step solution
1. Node 1 balance. \(T_0-2T_1+T_2=0\), or \(2T_1-T_2=120\).
2. Node 2 balance. \(T_1-2T_2+T_3=0\), or \(-T_1+2T_2=30\).
3. Solve. Simultaneous solution gives \(T_1=90^\circ\text{C}\) and \(T_2=60^\circ\text{C}\).
4. Residuals. \(r_1=120-2(90)+60=0\) and \(r_2=90-2(60)+30=0\).

Unit/sign/limit checks

Units: Every nodal equation is a temperature equation after division by \(kA/\Delta x\).

Signs: Temperatures decrease monotonically in the chosen positive heat-flow direction.

Limits: Both nodes lie between fixed boundary values and the zero residual confirms algebraic closure.

Final answer: \(\boxed{T_1=90^\circ\text{C},\quad T_2=60^\circ\text{C};\ r_1=r_2=0}\).

Engineering insight

The exact linear profile is a useful code-verification case: a conduction solver that cannot reproduce it has a mesh, indexing, or boundary-condition defect.

5–4

Two-Dimensional Steady Heat Conduction

On a square grid (\(\Delta x = \Delta y\)) with no generation, the steady energy balance gives an especially clean result: each interior node's temperature is the average of its four neighbors.

Interior node, 2-D square grid, no generation
\[ T_{\text{node}} = \frac{T_{\text{left}} + T_{\text{right}} + T_{\text{top}} + T_{\text{bottom}}}{4} \]
Interior node, 2-D square grid, with generation
\[ T_{\text{left}} + T_{\text{right}} + T_{\text{top}} + T_{\text{bottom}} - 4T_{\text{node}} + \frac{\dot{g}\,\Delta x^2}{k} = 0 \]

Solving the System

The resulting set of equations is solved directly (matrix inversion / Gaussian elimination) for small grids, or iteratively (Gauss–Seidel) for large ones, sweeping the grid repeatedly until temperatures stop changing.

Course Note

The "average of the four neighbors" rule is one of the most useful things to memorize. For a quick hand check of any 2-D steady problem with no generation, sweep the grid by hand using it a few times (Gauss–Seidel) — the numbers converge surprisingly fast.

EXAMPLE 2Interior Node in a Heat Spreader
Problem. A square-grid node in a ceramic heat spreader has neighboring temperatures \(82, 74, 91,\) and \(69^\circ\text{C}\). Uniform volumetric generation is \(6.0\times10^6\ \text{W/m}^3\), \(k=24\ \text{W/mK}\), and \(\Delta=4.0\ \text{mm}\). Determine the node temperature.

Strategy

An interior two-dimensional node with uniform spacing uses the four-neighbor balance plus the generation source.

Givens

  • Four neighbor temperatures; \(\dot g=6.0\times10^6\ \text{W/m}^3\).
  • \(k=24\ \text{W/mK}\), \(\Delta x=\Delta y=0.0040\ \text{m}\).

Physical/model interpretation

A unit-depth control volume conducts through four faces. Positive generation raises the node above the no-generation neighbor average.

Step-by-step solution
1. Formulate. \(\sum T_{\rm nbr}-4T_P+\dot g\Delta^2/k=0\).
2. Neighbor average. \((82+74+91+69)/4=79.0^\circ\text{C}\).
3. Source correction. \(\dot g\Delta^2/(4k)=6.0\times10^6(0.004)^2/[4(24)]=1.00\ \text{K}\).
4. Resolve. \(T_P=79.0+1.00=80.0^\circ\text{C}\).

Unit/sign/limit checks

Units: \((\text{W/m}^3)\text{m}^2/(\text{W/mK})=\text{K}\).

Signs: Positive generation produces a positive correction.

Limits: As \(\dot g\to0\), the answer returns to the four-neighbor average.

Final answer: \(\boxed{T_P=80.0^\circ\text{C}}\).

Engineering insight

The source correction scales with \(\Delta^2\), illustrating why coarse meshes can hide peak temperatures near concentrated electronic heat sources.

5–5

Transient Heat Conduction

For transient problems the time derivative is also discretized. Evaluating the spatial terms at the old time level gives the explicit method; at the new level, the implicit method.

Explicit scheme, interior node (1-D)
\[ T_m^{\,i+1} = \tau\left(T_{m-1}^{\,i} + T_{m+1}^{\,i}\right) + (1 - 2\tau)\,T_m^{\,i}, \qquad \tau = \frac{\alpha\,\Delta t}{\Delta x^2} \]

Stability Criterion (Explicit Method)

The explicit scheme is only stable when the coefficient of a node's own previous temperature stays non-negative. This caps the time step: in one dimension the mesh Fourier number must satisfy \(\tau \le \tfrac{1}{2}\) (and \(\tau \le \tfrac{1}{4}\) in two dimensions). Violating it produces physically impossible oscillating temperatures.

\[ \tau = \frac{\alpha\,\Delta t}{\Delta x^2} \le \frac{1}{2}\ \text{(1-D interior node)} \]

Explicit vs. Implicit

The explicit method computes each new nodal temperature directly but is limited by the stability criterion. The implicit method requires solving a system every step but is unconditionally stable, allowing larger time steps — the usual trade of cost per step against number of steps.

“Unconditionally stable” does not mean accurate: an implicit solution with an excessive time step may smear the transient. Explicit boundary nodes can impose a stricter criterion when convection or generation is present, so derive coefficients for the actual boundary equation.

EXAMPLE 3Stable Time Step for a Thermal-Shock Mesh
Problem. A one-dimensional glass slab model uses \(\Delta x=3.0\ \text{mm}\) and constant \(\alpha=5.5\times10^{-7}\ \text{m}^2/\text{s}\). An interior node and its neighbors initially read \(110,80,50^\circ\text{C}\). Find the largest stable explicit time step and the node temperature after one step using 80% of that limit.

Strategy

The requested direct time update identifies the explicit scheme; stability fixes the allowable mesh Fourier number before temperature is advanced.

Givens

  • \(\Delta x=0.0030\ \text{m}\), \(\alpha=5.5\times10^{-7}\ \text{m}^2/\text{s}\).
  • \(T_{m-1}^i=110^\circ\text{C}\), \(T_m^i=80^\circ\text{C}\), \(T_{m+1}^i=50^\circ\text{C}\).

Physical/model interpretation

Constant-property diffusion with no generation is modeled at an interior node. The symmetric neighbor temperatures imply zero discrete curvature.

Step-by-step solution
1. Stability limit. \(Fo_m=\alpha\Delta t/\Delta x^2\le0.5\), so \(\Delta t_{\max}=0.5\Delta x^2/\alpha\).
2. Calculate limit. \(\Delta t_{\max}=0.5(0.0030)^2/(5.5\times10^{-7})=8.18\ \text{s}\).
3. Choose step. \(\Delta t=0.80(8.18)=6.55\ \text{s}\), giving \(Fo_m=0.400\).
4. Advance. \(T_m^{i+1}=0.4(110+50)+(1-0.8)(80)=80.0^\circ\text{C}\).

Unit/sign/limit checks

Units: \(\Delta x^2/\alpha=\text{s}\); all update coefficients are dimensionless.

Signs: Coefficients \(0.4,0.2,0.4\) are nonnegative and sum to one.

Limits: The neighbor average equals \(80^\circ\text{C}\), so zero curvature correctly gives no local change.

Final answer: \(\boxed{\Delta t_{\max}=8.18\ \text{s};\ \Delta t=6.55\ \text{s},\ T_m^{i+1}=80.0^\circ\text{C}}\).

Engineering insight

Halving mesh spacing cuts the explicit stable step by four. Production models should demonstrate both spatial and temporal convergence rather than merely satisfying stability.

Chapter recap

  • Conservative control-volume balances expose the physics and handle boundaries cleanly.
  • Steady equations require residual and global-energy audits; transient equations additionally require stability and time-step checks.
  • Mesh refinement and credible properties are part of the solution, not optional presentation details.

Before moving on

  • Can you write interior, boundary, and source-node balances?
  • Can you distinguish stability, convergence, and accuracy?
← Chapter 4 · Transient Conduction Chapter 6 · Fundamentals of Convection →