# Chapter 2 One-Dimensional Steady State Conduction Analytical and Numerical Methods

## 1  One-Dimensional Conduction in Cartesian Coordinates

### 1.1  Plane Wall, Dirichlet and Newmann Conditions

Consider a flat wall of homogeneous, isotropic material, with thermal conductivity k extending from x=X1 to x=X2.
The heat equation is
 - Ñ·q + g = - dq dx + g = 0
where
 g = d (Q/V) dt = dQV dt = × QV
is the rate of internal thermal energy generation per unit volume.
For constant conductivity and no internal heat generation (g=0)
 d2 T dx2 = 0
Boundary conditions (Dirichlet):
 T(x=X1) = T1
 T(x=X2) = T2
Boundary conditions (Dirichlet-Newmann):
 T(x=X1) = T1
 -k d T dx |x=X2 = q2
Boundary conditions (Dirichlet-Robin):
 T(x=X1) = T1
 -k d T dx |x=X2 = h2 (T2 - T¥)
 T(x=X1) = T1
 -k d T dx |x=X2 = es(T24 - T¥4)

## 2  One-Dimensional Conduction in Polar Coordinates

### 2.1  Cylindrical Shell

Consider a hollow cylindrical shell of homogeneous, isotropic material, with thermal conductivity k extending from r=R1 to r=R2.
Heat equation:
 - Ñ·q + g = - 1 r d dr (r q) + g = 0
For constant conductivity and no internal heat generation (g=0)
 d dr (r d T dr ) = d2 T dr2 + 1 r dT dr = 0
Boundary conditions (Dirichlet):
 T(r=R1) = T1
 T(r=R2) = T2
Boundary conditions (Dirichlet-Newmann):
 T(r=R1) = T1
 -k d T dr |r=R2 = q2
Boundary conditions (Dirichlet-Robin):
 T(r=R1) = T1
 -k d T dr |r=R2 = h2 (T2 - T¥)
 T(r=R1) = T1
 -k d T dr |r=R2 = es(T24 - T¥4)

### 2.2  Spherical Shell

Hollow spherical shell, homogeneous, isotropic material, thermal conductivity k extending from r=R1 to r=R2.
Heat equation:
 - Ñ·q = g = - 1 r2 d dr (r2 q) + g = 0
For constant conductivity and no internal heat generation (g=0)
 d dr (r2 d T dr ) = d2 T dr2 + 2 r dT dr = 0
Boundary conditions (Dirichlet):
 T(r=R1) = T1
 T(r=R2) = T2
Boundary conditions (Dirichlet-Newmann):
 T(r=R1) = T2
 -k d T dr |r=R2 = q2
Boundary conditions (Dirichlet-Robin):
 T(r=R1) = T1
 -k d T dr |r=R2 = h2 (T2 - T¥)
 T(r=R1) = T1
 -k d T dr |r=R2 = es(T24 - T¥4)

## 3  Concept of Thermal Resistance

The heat flux q and rate of heat flow [Q\dot] are related by
q = d (Q/A)

dt
=
×
QA

=
 × Q

A
= -k ÑT
Therefore
 × Q = -k A ÑT
In Cartesian coordinates for a wall of span X2-X1 = L subject to Dirichlet conditions such that T2 - T1 = DT, such that the flow of heat is along the positive x- direction,
 × Q = -k A dT dx = -k A DT L = - DT R
where the thermal resistance R is defined as
 R = L A k
In Cylindrical coordinates for a shell of span R2 - R1 = L, subject to Dirichlet conditions such that T2 - T1 = DT, such that the flow of heat is along the positive r- direction
 × Q = - DT R
where the thermal resistance R is defined as
 R = ln(R2/R1) 2 pk L
In Spherical coordinates for a shell of span R2 - R1 = L, subject to Dirichlet conditions such that T2 - T1 = DT, such that the flow of heat is along the positive r- direction
 × Q = - DT R
where the thermal resistance R is
 R = R2-R1 4 pk ra rb
Since at a solid-fluid interface q = [Q\dot]/A = h (Ts - T¥) the convective thermal resistance there is
 R = 1 h A
Thermal resistances are analogous to electrical resistances and for heat flowing through a series of resistances, the overall thermal resistance is obtained by simple addition. To illustrate, consider a plane wall subjected to convective heat exchange at surfaces X1 and X2, with heat exchange coefficients and fluid temperatures, respectively given by h1, T1,¥, h2, T2,¥, the rate of heat flow is given by
 × Q = - DT Ro = - Tb,¥ - Ta,¥ Ra + RL + Rb
where Ra = 1/ha A, RL = L/k A and Rb = 1/hb A.
Sometimes, the concept of overall heat transfer coefficient U is used. It is related to the overall thermal resistance by
 U = 1 A Ro
An important related concept which will be useful later on is the Biot number Bi. It is a measure of the relative magnitude of internal and external resistances to heat flow and is defined as
 Bi = h L k
The above ideas are useful when studying steady state conduction through composite walls or shells. A direct application is the notion of critical radius of insulation. Consider a long pipe carrying a hot fluid and coated with insulation of thermal conductivity k. Taking the derivative of the overall thermal resistance with respect to the outer radius of insulation for constant h at the outer radius of insulation yields
 d Ro dr = 1 2 pr ( 1 k - 1 hr )
Note the overall thermal resistance has an extremum at
 rc = k h
Taking the second derivative one can show the extremum is minimum.

## 4  Problems Involving Internal Heat Generation

For problems involving non-zero internal heat generation, the governing differential equations are simply obtained by maintaining the term g in the equation.
Specifically, for Cartesian coordinates
 - d dx (k dT dx ) = g
For cylindrical coordinates
 - 1 r d dr (r k dT dr ) = g
And for spherical coordinates
 - 1 r2 d dr (r2 k dT dr ) = g

## 5  Numerical Methods

Often, analytical solutions are not possible and one must use numerical approximation techniques. The most commonly used numerical approximation methods are: finite differences, finite volumes and finite element methods. They are all based on the conversion of the original problem into a discrete, algebraic form. Standard techniques of numerical algebra are then used to obtain the desired approximation.

### 5.1  Finite Difference/Finite Volume Method

In the finite difference method one starts with the differential formulation and by a process of discretization transforms the problem into a system of interlinked simultaneous algebraic equations that then must be solved in order to determine an approximation to the desired solution. Discretization consists of first introducing a mesh of nodes by subdividing the solution domain into a finite number of sub domains and then approximating the derivatives in the boundary value problem by means of appropriate finite difference ratios which can be obtained from a truncated Taylor series expansion Finite difference method) or by integration over small volume elements (finite volume method). A system of interlinked simultaneous algebraic equations is obtained that then must be solved in order to determine an approximation to the desired solution.
Consider the following steady state heat conduction problem inside a homogeneous, isotropic flat wall of span L = X2-X1, constant thermal conductivity, with internal heat generation
 -k d2 T dx2 = g
and subjected to Dirichlet conditions
 T(x=X1) = 0
 T(x=X2) = 0
Introduce a mesh in [X1,X2] by dividing the interval into N equal subintervals of size h. This produces two boundary mesh points x1=X1 and xN+1=X2, and N-1 interior mesh points xi, i=2,..,N-1 for a total of N+1 mesh points. Since the values T(x1) = 0 and T(xN+1) = 0 are known, the values T(xi), i=2,..,N are to be determined.
From the Taylor expansions for T(xi+1) and T(xi-1) the following centered difference formula is obtained
 d2 T dx2 |xi = 1 h2 [T(xi+1) - 2 T(xi) + T(xi-1)] + h2 12 d4 T dx4 (xi)
for some xi Î (xi-1,xi+1).
If the higher order terms are discarded from the above formula and the approximations are introduced in the original differential equation, the result, after slight rearrangement and with Ti = T(xi), is
 - Ti-1 + 2 Ti - Ti+1 = h2 ( g k )
This is a tridiagonal system of simultaneous linear algebraic equations for the unknown approximate values of Txi = Ti, for i=2,3,...,N. In matrix notation the system becomes
 A u = b
where A is a tridiagonal matrix, a sparse, positive definite, diagonally dominant matrix with very nice properties and, for the specific boundary value problem considered above, given (in augmented form) by
A = æ
ç
ç
ç
ç
ç
ç
ç
ç
ç
ç
è
 1
 0
 0
 0
 0
 ...
 0
 -1
 2
 -1
 0
 0
 ...
 0
 0
 -1
 2
 -1
 0
 ...
 0
 ..
 ..
 ..
 ..
 ..
 ...
 0
 ..
 ..
 ..
 ..
 ..
 ...
 0
 0
 0
 0
 -1
 2
 -1
 0
 0
 0
 0
 0
 -1
 2
 -1
 0
 0
 0
 0
 0
 ...
 1
ö
÷
÷
÷
÷
÷
÷
÷
÷
÷
÷
ø
Moreover, u is the (column) vector of unknowns given in this case by
u = æ
ç
ç
ç
ç
ç
ç
ç
ç
ç
è
 T1
 T2
 T3
 ...
 ...
 TN
 TN+1
ö
÷
÷
÷
÷
÷
÷
÷
÷
÷
ø
and b is the (column) forcing vecto, again, given in this case by
f = h2 ( g

k
) æ
ç
ç
ç
ç
ç
ç
ç
ç
ç
è
 0
 1
 1
 ...
 ...
 1
 0
ö
÷
÷
÷
÷
÷
÷
÷
÷
÷
ø
An identical procedure is used to obtain the corresponding algebraic set of equations for problems in polar coordinates.

### 5.2  Finite Element Method

In the finite element method, one starts with the variational formulation and then, by a process of finite element function representation transforms the problem into a system of interlinked simultaneous algebraic equations.This system must be solved in order to determine an approximation to the desired solution.
Consider again the problem of steady state heat conduction through a flat wall of span L = X2-X1 = 1 - 0 = 1 and of uniform, isotropic material of constant thermal conductivity, with internal energy generation and subject to Dirichlet boundary conditions. To keep things simple, assume the boundary conditions are homogeneous, i.e.
 T(x=0) = 0
 T(x=1) = 0
The variational formulation is obtained by introducing a function v(x) satisfying v(0)=0 and ò01 v2 dx < ¥, but otherwise rather arbitrary and then integrating the differential equation, i.e.
 óõ 1 0 - d2 T dx2 v dx = óõ 1 0 g k v dx
Integrating the left hand side by parts readily yields the variational or weak form
 a(T,v) = (f,v)
where
 a(T,v) = óõ 1 0 dT dx dv dx dx
and
 (f,v) = óõ 1 0 f v dx
where f = g/k.
In the Galerkin finite element method, one seeks an approximate solution directly using the variational statement of the problem a(T,v) = (f,v). One starts by introducing a set of basis functions {fi(x)}i=1N satisfying the stated boundary conditions. The number N is the number of finite elements in the system.
The approximate solution is then expressed as a simple linear combination of the basis functions by
 TN(x) = Nå i=1 Ti fi(x)
Note that when using this representation, the quantities Ti are simply the numerical coefficients of the basis functions in the representation and are independent of the variable x. Most importantly, these coefficients turn out to be the values of the computed approximation at the nodal locations of the finite element mesh. Furthermore, the test function v(x) in this case is selected as a simple linear combination of the basis functions,
 v(x) = Nå j=1 fj(x)
The above choices for TN(x) and v(x) are introduced into the variational statement in order to obtain the discrete system.
The finite elements are introduced when the given domain x Î [X1,X2] is subdivided into a collection of N contiguous subdomains [x1=X1,x2], [x2,x3],...,[xi,xi+1], ..., [xN,xN+1=X2]. The points x1, x2,x3,...,xi,...,xN+1 are called nodes and each contiguous pair of them, [xi,xi+1], represents the boundaries of the element that spans the space hi = xi+1-xi between the nodes (nodal spacing). In this case, each element has two nodes one located at each end of the element. The location of all these points is given with reference to a global system of coordinates. If nodes are numbered based on their position is terms of global coordinates one obtains global node numbers. For simplicity, assume all finite elements are of equal size (uniformly spaced mesh of nodes) so that h1 = h2 = ... = hi = h.
In the local description of the element bounded by nodes [xi,xi+1], the positions measured with respect to the global system of coordinates are transformed into coordinates [x1,x2] = [-1,+1] referred to a system of coordinates with its origin in the center of the element, i.e. The affine or linear transformation is simplest, i.e.
 x(x) = 2 x - xi - xi+1 hi
This is called the standard element. The inverse transformation is
 x(x) = 1-x 2 xi + 1+x 2 xi+1 = y1(x) xi + y2(x) xi+1
Note that the functions y1(x) and y2(x) can be readily written as functions of the global coordinate x by substituting the affine mapping for x, i.e.
 y1(x) = 1 2 ( 1- 2 x - xi - xi+1 hi )
 y2(x) = 1 2 ( 1+ 2 x - xi - xi+1 hi )
Global basis functions fi(x) (i.e. basis functions in the global coordinate system) for each nodal location xi where elements ei-1 and ei meet, are readily obtained from the above expressions for the standard element as follows
fi(x) = ì
ï
í
ï
î
 y2(x)
 if xi-1 £ x £ xi
 y1(x)
 if xi £ x £ xi+1
 0
 otherwise
The role of the local basis functions is to generate the value of the approximated quantity inside the element from the values at the nodes by interpolation. For instance, for an arbitrary interior element ei where the index i denotes the global node number. If the values of the required solution at xi and xi+1 are u(xi) = ui and u(xi+1) = ui+1, respectively, the approximate values of T(x), T(x)ei inside the element are calculated by simple interpolation as follows
 Tei(x) = Ti f1i + Ti+1 f2i
Here, the local finite element basis functions for element ei, f1i and f2i are defined as
 f1i = xi+1 - x hi

 f2i = x - xi hi
where hi = xi+1 - xi is the nodal spacing (and, in this case, also the element size) and all the positions are measured in the global coordinate system. Note that the indices 1 and 2 on the local basis functions refer to the local node numbers for the element.
Associated with each node of the finite element mesh, there is a global basis function fi(x), i=1,2,...N which is readily constructed once the local finite element basis functions are specified.
The approximate solution on the entire domain Th(x) is represented as a linear combination of the global finite element basis functions, i.e.
 Th(x) = T1 f1(x) + T2 f2(x) + ... + Ti fi(x) + TN+1 fN+1(x) = N+1å i=1 Ti fi(x)
The specific relationships between local and global finite element basis functions are obtained as follows. For node 1,
f1 = ì
ï
í
ï
î
 f11 = x2-x x2-x1
 if x Î e1 = [x1,x2]
 0            elsewhere
for node i,
fi = ì
ï
ï
ï
í
ï
ï
ï
î
 f2i-1 = x-xi-1 xi-xi-1
 if x Î ei-1=[xi-1,xi]
 f1i = xi+1-x xi+1-xi
 if x Î ei=[xi,xi+1]
 0 elsewhere
and for node N+1
fN+1 = ì
ï
í
ï
î
 f2N+1 = x-xN xN+1-xN
 if x Î eN=[xN,xN+1]
 0            elsewhere
Finally, as mentioned before, in the Galerkin formulation of the finite element method, the test functions v(x) are selected as
 v(x) = N+1å i=1 fi(x)
Introducing all the above into the variational formulation of the problem and rearranging (with, for simplicity X1=0 and X2=1, yields the Galerkin finite element method equations as
 Tj-1 óõ xj+1 xj-1 d fj-1 dx d fj dx dx + Tj óõ xj+1 xj-1 ( d fj dx )2 dx + Tj+1 óõ xj+1 xj-1 d fj+1 dx d fj dx dx = óõ xj+1 xj-1 f fj dx
for all j=1,3...N+1, or more succinctly as
 Nå i=1 Ti a(fi, fj) = (f, fj)
for all i,j=1,2,...,N+1. Here
 a(fi, fj) = óõ 1 0 d fi dx d fj dx dx
and
 (f, fj) = óõ 1 0 f fj dx
for i,j=1,2,...,N+1.
Using matrix notation, the above is simply written as
 K u = F
Solving the algebraic problem yields the values of the nodal values ui, i=1,2,...,N+1. The matrix K is called the finite element stiffness matrix, the colum vector F is the force vector and the column vector u is the vector of unknown nodal values of u.

## 6  Extended Surfaces

Extended surfaces are widely used to enhance heat exchange. By attaching fins to a hot wall, heat is more easily dissipated and the wall temperature is reduced. The enhanced heat dissipation takes place through the extended surface of the fin. Consider a thin fin attached to a hot flat wall at x=0. The fin has cross sectional area A(x) and perimeter P(x). Since the fin is thin, only temperature variations along the length of the fin will be considered. Consider a differential element in the fin of thickness Dx located at some distance x from the base of the fin. Under steady state conditions, of the heat entering the element a fraction is dissipated by convection through the extended surface and the rest is transferred further down along the fin. The result of writing the down the differential thermal energy balance and computing the limit as Dx ® 0 is
 d dx (A(x) dT dx ) - h Px k (T - T¥) = 0
Once A(x), P(x), h, k, and T¥ are specified, together with suitable boundary conditions, the differential equation can be solved to obtain the temperature distribution along the fin T(x).
The heat transfer enhancement effectiveness of a fin is assessed by computing its fin efficiency. This is defined as the ratio of the actual amount of heat convected through the extended surface of the fin to the convected heat that would obtain if the fin were at the base temperature throughout.

File translated from TEX by TTH, version 3.85.
On 23 Sep 2009, 17:13.

Updated: 2009-09-23, 17:13