Introduction to Numerical Methods

in Heat Transfer

Presented by:
Steven L. Rickman
NASA Technical Fellow for Passive Thermal
NASA Engineering and Safety Center

Thermal and Fluids Analysis Workshop

Cleveland, Ohio
August 2024
1
NASA Photo: ISS016E009207
 Introduction

Heat transfer is best understood through theory and

application of principles in thermal analysis;

Modern thermal analysis leverages the power of

computers and numerical methods to simulate heat
transfer in networks representing a physical system;

This lesson is an introduction to numerical methods

in heat transfer.

2
 Overview

Students may have experience with numerical

methods in college level courses;

Once in the workplace, however, they often use

commercially available modeling tools;

Unless graduate study is pursued, students must

create opportunities to understand the
fundamentals behind these methods;

This course is intended to provide such focus.

3
 Scope of this Lesson

What is heat transfer?

Governing differential equation;

Finite difference;

Finite element;

Thermal radiation in heat transfer analysis.

4
 Lesson Roadmap

Thermal Fourier's Heat

What is Heat Transfer?
Modeling Basics Conduction Law

Governing Differential Example 1: Solving the

Solution for 1-D Heat
Equation for 1-D Heat Transient Heat Equation
Transfer
Transfer (1-D)

Example 2: Example 3:
Finite Implicit
Transient Thermal Adding Thermal
Difference Formulation
Numerical Analysis Radiation
Methods
Establishing the
Finite Calculus of
Functional
Element Variations
Equation

Example 5:
Example 4: Steady State Analysis Example 6: Transient
Steady State
- Distributed Heating Analysis
Analysis

5
 What is Heat Transfer? (Ref. 1)

Heat transfer is energy transfer due to a temperature

difference and can only be measured at the boundary
of a system.

Conduction - Heat transfer from one substance to

another by direct contact.

Convection - Heat transfer via movement of fluids.

Radiation – Heat transfer via electromagnetic waves.

6
 Thermal Modeling Basics

Engineers use thermal models as a tool to aid in

design and to understand the performance of
thermal systems;

A thermal model is an abstraction of a physical

system.

7
 Thermal Modeling Basics

Closed-form solutions do not exist for all physical

systems or geometries -- so they are discretized
into smaller pieces, called elements.

Geometry Geometry Discretized into a

Note: Geometry and mesh created using MSC Patran®
Finite Element Mesh
8
 Thermal Modeling Basics

Thermal systems can also be represented by

representations analogous to electrical resistance-
capacitance (RC) networks.

Geometry that is discretized is represented by

nodes. A node is isothermal and has a constant
temperature throughout the entire volume.

Heat flow between nodes is modeled with

conductors (the inverse of thermal resistance) and
heat storage is modeled using a capacitor.
9
 Thermal Modeling Basics

Nodes come in the following varieties:

Diffusion -- a node representing a finite mass with

finite capacitance;

Arithmetic -- a massless node representing zero

capacitance;

Boundary -- a node representing a source or sink

with infinite capacitance.

10
 Thermal Modeling Basics

Conductors come in the following varieties:

Conduction -- heat transfer between solid objects

(or a solid and a gas/fluid);

Convection -- heat transfer between a solid object

and a convecting liquid or gas;

Radiation -- heat transfer via electromagnetic

radiation between objects.

11
 Thermal Modeling Basics

Heat transfer via conduction and convection is

proportional to DT:

But, heat transfer via radiation takes this form:

12
 Thermal Modeling Basics

Consider this
geometry

The network representation for this configuration

might look something like this:

13
 Some Terminology

Heat rate, denoted by , is [energy/time];

Heat flux, denoted by , is [energy/area/time];

Volumetric heat rate, denoted by , is

[energy/volume/time];

Heat per unit length, denoted by , is

[energy/length/time].

14
 Fourier Heat Conduction Law

Fourier assumed that

conduction heat
transfer is directly T1 T2
proportional to the k
A A
temperature difference
across a solid:
L

15
 Conservation of Energy

In heat transfer calculations, energy must be

conserved:

where...
is energy per unit time entering;
is energy per unit time generated;
is energy per unit time leaving;
is time rate of change of energy stored.
16
Deriving the Heat Equation in One Dimension
(Ref. 2)

x x+dx

17
 Deriving the Heat Equation in One Dimension
(Ref. 2)
The energy per unit time stored in a mass can be
expressed as:

where...
is the mass [mass];
is the specific heat [energy/mass/temperature];
is the material density [mass/volume];
is the material volume [volume];
is the differential element cross sectional area [area];
is the differential element length [length];
is the mass' time rate change of temperature [temperature/time].
18
 Deriving the Heat Equation in One Dimension
(Ref. 2)
The heat entering the mass per unit time at location x
is:

where...
is the thermal conductivity [energy/time/length/temperature];

is the differential element cross sectional area [area];

is temperature gradient at x [temperature/length].

19
 Deriving the Heat Equation in One Dimension
(Ref. 2)
The heat leaving the mass per unit time at location
x+Dx is:

where...
is the thermal conductivity [energy/time/length/temperature];

is the differential element cross sectional area [area];

is temperature gradient at x+dx [temperature/length].

20
 Deriving the Heat Equation in One Dimension
(Ref. 2)
But we recognize that the heat transfer at x+dx may
be expressed as:

where, the first term is simply the heat entering the

left face per unit time and the second term is the
change in heat transfer rate with respect to x over the
control volume times the distance to the other end of
the control volume.
21
 Deriving the Heat Equation in One Dimension
(Ref. 2)
We can now substitute our new expressions for the
original terms:

This reduces to:

22
 Deriving the Heat Equation in One Dimension
(Ref. 2)
Rearranging and simplifying yields:

Or, for a constant thermal conductivity:

This is known as the heat equation in one

dimension expressed per unit volume.
23
 Extension to 2-D and 3-D

More generally, the heat equation for isotropic

thermal conductivity can be expanded into
additional dimensions:

or, using the Laplacian:

24
 One Dimensional Heat Equation:
Steady State
Assuming the heat generation term is constant, the
only time dependency appears on the right hand
side of the equation: 0

At steady state, the time rate of change of

temperature, ∂T/ ∂ t = 0;

We see that steady state behavior is independent of

density, r, and specific heat, Cp .
25
 One Dimensional Heat Equation:
Transient
For a case with no internal heat generation, we
have:
0

We recognize this to be a linear differential equation

of second order in the spatial dimension, x and first
order in time, t.

Fortunately, there exists a solution.

26
 Example 1: Solving the Transient Heat
Equation (Adapted from Ref. 2)
Consider the following example:
r , Cp , k

T1 T0 T1
x=0 x=2L
x
The rod is initially at a uniform temperature, T0. At
time, t = 0, the temperature at both ends (x=0 and
x=2L) is raised to T1;

How does the temperature profile in the rod change

with time? 27
 Example 1: Solving the Transient Heat
Equation (Adapted from Ref. 2)
The governing differential equation is the previously
derived heat equation:

We recast the equation noting that a = k/rCp and

q=T-T1:

28
 Example 1: Solving the Transient Heat
Equation (Adapted from Ref. 2)
Using this new form, we examine our boundary
conditions for q(x,t):

For 0 < x < 2L and t = 0: q(x,0) = T0 - T1

For x = 0 and t > 0: q(0,t) = 0

For x = 2L and t > 0: q(2L,t) = 0

29
 Example 1: Solving the Transient Heat
Equation (Adapted from Ref. 2)
To solve this differential equation, we can use
separation of variables by assuming:

where...

X(x) is a function of only x, and;

Y(t) is a function of only t.

30
 Example 1: Solving the Transient Heat
Equation (Adapted from Ref. 2)
Forming the required derivatives of q(x, t):

So our transformed differential equation becomes:

31
 Example 1: Solving the Transient Heat
Equation (Adapted from Ref. 2)
We can rearrange the previous equation to read:

To facilitate solution, we can set both equations equal

to a constant, where l is called the separation
constant:

32
 Example 1: Solving the Transient Heat
Equation (Adapted from Ref. 2)
We can now write two, independent linear
differential equations:

33
 Example 1: Solving the Transient Heat
Equation (Adapted from Ref. 2)
We recognize that the first equation:

has a solution of the form:

34
 Example 1: Solving the Transient Heat
Equation (Adapted from Ref. 2)
The second equation:

has a solution of the form:

35
 Example 1: Solving the Transient Heat
Equation (Adapted from Ref. 2)
Our expression for q(x, t), then, becomes:

We can simplify this expression by letting:

36
 Example 1: Solving the Transient Heat
Equation (Adapted from Ref. 2)
The expression becomes:

The previously defined boundary conditions are

now used to solve for constants, C1 and C2:

For 0 < x < 2L and t = 0: q(x,0) = T0 - T1

For x = 0 and t > 0: q(0,t) = 0
For x = 2L and t > 0: q(2L,t) = 0
37
 Example 1: Solving the Transient Heat
Equation (Adapted from Ref. 2)
When we apply the second boundary condition (i.e.,
for x = 0 and t > 0):

From this, we see that C1 must be zero.

We note that C2 cannot be zero -- if it was, q(x, t) = 0

everywhere.
38
 Example 1: Solving the Transient Heat
Equation (Adapted from Ref. 2)
When we apply the third boundary condition (i.e., x
= 2L and t > 0):

Since C2 ≠ 0, we see that the only way for this to

happen is when:

39
 Example 1: Solving the Transient Heat
Equation (Adapted from Ref. 2)
This happens when:

So the solution may be expressed in the form of a

series given by:

40
 Example 1: Solving the Transient Heat
Equation (Adapted from Ref. 2)
Cn may be determined by integrating considering
the initial conditions (for n = 1, 3, 5, ...):

The overall solution becomes (for n = 1, 3, 5, ...):

41
 Example 1: Solving the Transient Heat
Equation
Time = 0.1 seconds
T0 = 100 C
T1 = 200 C
n=1
n=3
Temperature (C)

n=5
n = 21

Fraction of Rod Length (x/2L)

42
 Example 1: Solving the Transient Heat
Equation
n = 199
T0 = 100 C
T1 = 200 C
Temperature (C)

t = 1 sec
t = 2 sec
t = 5 sec
t = 10 sec
t = 20 sec
t = 100 sec

Fraction of Rod Length (x/2L)

43
 Numerical Methods

We'll focus on two different numerical methods:

Finite Difference -- uses the differential formulation

-- i.e., equations are formulated using the governing
differential equation -- where we replace the partial
derivatives by approximations obtained by Taylor
expansions near the point of interest;

Finite Element -- uses a variational formulation --

i.e., equations are formulated from an integral
formulation arising from the Calculus of Variations.
44
 Formulation of the Finite Difference for 1-D
Heat Transfer
Finite difference relies on a differential formulation -
- that is, a description of the heat transfer using
derivatives;

For our one-dimensional heat transfer case, recall

the governing differential equation is:

45
 Formulation of the Finite Difference for 1-D
Heat Transfer
This equation involves the second derivative of
temperature with respect to a spatial dimension
(i.e., x) and the first derivative of temperature with
respect to time;

Let's take a closer look at how these derivatives are

formed for a numerical solution.

46
 Formulation of the Finite Difference for 1-D
Heat Transfer
Consider this temperature distribution about x:

Dx Dx
T(x)

x -Dx x x +Dx

x
47
 Formulation of the Finite Difference for 1-D
Heat Transfer (Ref. 3)
At a given instant in time, the first derivative, taken
to the "right" of x, is given by:

Similarly, the first derivative, taken to the "left" of x,

is given by:

Each expression yields an approximation of the

slope in the specified region about x.
48
 Formulation of the Finite Difference for 1-D
Heat Transfer (Ref. 3)
The second derivative is just the slope of the slope
over the region, or, the slope of the first derivatives:

which becomes...

49
 Formulation of the Finite Difference for 1-D
Heat Transfer (Ref. 3)
Further simplification and rearranging yields:

Hence we have an expression for the second

derivative in terms of temperatures at specific
locations.

But this derivation assumed the distance between

nodes was Dx on either side of the node. What if it
isn't?
50
 Formulation of the Finite Difference for 1-D
Heat Transfer (Adapted from Ref. 3)
Consider the case where the distance between the
last node and the boundary is only Dx/2:

If we go through the same process, our expression

for the second derivative at the left end and right
end, respectively becomes:

51
 Formulation of the Finite Difference for 1-D
Heat Transfer (Ref. 3)
Forming the first derivative of temperature, T with
respect to time, t for a specific node is considerably
easier:

Again, we have an expression for a derivative in

terms of temperature and time.

52
 Subscripting and Superscripting

We're going to be working with both time and

spatial dimension;

We'll need to do some bookkeeping;

A subscripting and superscripting scheme is shown

here and will be used in the example.

53
 Subscripting and Superscripting

For the node of interest, i, at time, n, we have:

i-1 i i+1

Time T of T of T of
Node i-1 Node i Node i+1
n
n+1

54
 Review

We have...

derived expressions for the second derivative of

temperature with respect to distance and for the
derivative of temperature with respect to time, and;

established a subscript and superscript convention

to aid our analysis;

We're now ready to present an example problem.

55
 Example 2: Explicit Finite Difference

Consider the following configuration:

20 C 100 W
(for t > 0)
20 C
(const.) aluminum (const.)

r = 0.01 m
L = 0.1 m

For an initial rod temperature of 20 C, determine

the system transient response.
56
 Example 2: Explicit Finite Difference

We'll model the system with five diffusion ( ) nodes:

0.01 m 100 W 0.02 m

20 C (const.) 20 C (const.)
1 2 3 4 5

L = 0.1 m

We also establish two boundary ( ) nodes at either

end of the rod.
57
 Example 2: The Explicit Finite Difference
Solution
Consider the following simplified geometry;

We wish to write the difference equation for node i

accounting for heat transfer via conduction to node
i+1 and the addition of heat to node i.

58
 Example 2: The Explicit Finite Difference
Solution
The overall energy balance for a segment of the rod
is:

where each term represents energy per unit time

per unit volume;

This leads to the governing differential equation:

59
 Example 2: The Explicit Finite Difference
Solution
We arrive at the difference equation describing our
system:

which becomes...

For an explicit scheme, note that the temperature at

time, n+1 is completely determined by parameters
known at time n. 60
 Example 2: The Explicit Finite Difference
Solution
Temperatures at time n are known and are used to
calculate temperatures at time n+1 -- the n+1
solution becomes the n solution for the next
iteration.

61
 Example 2: Explicit Finite Difference
70
x= 0.00 m x = 0.01 m
x= 0.03 m x = 0.05 m
x= 0.07 m x = 0.09 m
60 x= 0.10 m
Temperature (C)

50

40

30

20
0 5 10 15 20 25 30 35 40
Time (s)
62
 Example 2: Explicit Finite Difference
70
t=0 s
60 t=5s
t = 10 s
Temperature (C)

t = 20 s
50
t = 150 s

40

30

20

10
0.00 0.02 0.04 0.06 0.08 0.10

X Location (m) 63
 Time Step and Time Constant

Finite difference temperature solution is calculated

at discrete time steps;

Accuracy and stability of the solution is influenced

by the time step;

Commercial software will auto calculate the time

step;

Upper limit to the time step required for solution

stability- function of the system's time constant, t.
64
 Time Step and Time Constant

For transient solutions, the time constant, t, is

determined by:

In other words, the minimum of the quotient of

node capacitance divided by the sum of the
conductances to that node is the limiting rate at
which marching in time can occur.

65
 Example 2: Effect of the Time Step

Time Step = 1.0 s

Time Step = 2.5 s

Time Step, t < 1.93 s,

required for stability
Time Step = 3.5 s
66
 Adding Radiation

Conduction and convection heat transfer are linear

functions of the temperature difference;

Radiation heat transfer is nonlinear and is

proportional to difference of the fourth power of
absolute temperatures.

67
 Adding Radiation (Ref. 5)

We seek to express the radiation heat transfer

between the two nodes of interest in terms of DT.

Linearized Radiation Conductor

68
 Adding Radiation

Consider the following simplified geometry;

We wish to write the difference equation for node i

including radiation to the environment.

69
 Adding Radiation

We can perform an energy balance on node i

(adapted from Ref. 6):

This is an expression of energy transfer per unit time

per unit volume;

70
 Adding Radiation

Expanding and assuming nodes are of equal length,

Dx:

Add nodal subscripting and linearization:

71
 Adding Radiation

Next, we add accommodation for node-specific

heating and a means of indexing the equation in
time (i.e., n, n+1):

With some additional manipulation, we arrive at our

desired result:

72
 Example 3: Explicit Finite Difference with
Radiation Added
Consider the following configuration:
Tenv = 20 C, -250 C
20 C 100 W e= 0.8 20 C
(const.) aluminum (const.)

r = 0.01 m
L = 0.1 m

Assume a surface e = 0.8 and solve for two different

Tenv: 20 C and -250 C.
73
 Example 3: Explicit Finite Difference with
Radiation Added
70
Steady State Temperature (C)

60

50

40
No Radiation

30 With Radiation (20 C Boundary)

With Radiation (-250 C Boundary)

20
0.00 0.02 0.04 0.06 0.08 0.10
X Location (m)
74
 The Implicit Finite Difference Solution

The explicit technique shown allows determination

of temperatures at time n+1 in terms of known
temperatures at time n;

However, different formulations of the differencing

equation exist;

Common Implicit formulations include the Backward

difference and the Central difference (Crank-
Nicolson) methods.

75
 The Implicit Finite Difference Solution
(Adapted from Ref. 3)
Crank-Nicolson takes the form:

Note that since the n+1 superscript appears on both

sides of the equation, we can't express
temperatures at that time as explicit functions of n;

Hence, this is referred to as an implicit technique.

76
 Solution Accuracy (Adapted from Ref. 4)

An assessment of the solution accuracy can be

made by observing the order of the terms truncated
in the Taylor series approximation;

For a first derivative formed as:

The truncated terms are:

77
 Solution Accuracy (Adapted from Ref. 4)

A second derivative with respect to x, though, is

inherently first order accurate in x because it is
formed using first derivatives that are accurate to
the first order in x.

Techniques exist to boost Forward and Backward

differences to higher accuracy.

78
 Solution Accuracy (Adapted from Ref. 4)

Both Forward and Backward differencing are

inherently first order accurate in x;

A Central difference, such as Crank-Nicolson,,

however, is inherently:

79
 Finite Difference Wrap-Up

Formulated the finite difference;

Demonstrated explicit solution technique;

Brief discussion of time step and time constant;

Added effects of radiation;

Implicit techniques;

Solution accuracy.
80
 Finite Element Strategy

The finite element formulation centers around

minimization of a functional;

But what is a functional?

To understand this, we need some background on

the Calculus of Variations.

81
 Calculus of Variations

Finite differencing relies on a differential

formulation of the heat equation;

Finite element relies on a variational formulation;

For steady state, one-dimensional heat transfer, we

seek to minimize the following integral, called a
functional:

82
 Calculus of Variations (Adapted from Ref. 7)

But where does the functional come from and what

is the theoretical basis for establishing the
functional?

The functional is established using the Calculus of

Variations;

In general, we seek a function to minimize the

integral:

83
 Calculus of Variations (Adapted from Ref. 7)

u(x) is the function we want to minimize.

ũ(x,e) is the set of all functions that can minimize
the integral, I.

ũ(x,e)

e(x)
u(x)

x=0 x=L
84
 Calculus of Variations (Adapted from Ref. 7)

The set of all functions to minimize I can be

expressed as:

where...
u(x) is the function we seek to minimize;
(x) is an arbitrary function constrained by:

85
 Calculus of Variations (Adapted from Ref. 7)

This ensures that ũ(x,e) is correct at the boundaries:

This works if e = 0 so that:

or... 0

86
 Calculus of Variations (Adapted from Ref. 7)

Let's substitute our expression for ũ(x,e) into the

original equation but note we've added the
variable, e :

and note that when e = 0, we return to the original

expression.

87
 Calculus of Variations (Adapted from Ref. 7)

So, we want our function I(e) to be a minimum

when e = 0;

We can accomplish this by differentiating the

expression with respect to e;

88
 Calculus of Variations (Adapted from Ref. 7)

By chain rule, the expression becomes:

This looks messy, but we can clean it up.

89
 Calculus of Variations (Adapted from Ref. 7)

Remember:

So that...

and...

90
 Calculus of Variations (Adapted from Ref. 7)

The integral expression becomes:

But note that:

w dv
can be integrated by parts.
91
 Calculus of Variations (Adapted from Ref. 7)

Recall integration by parts:

So the entire integral expression becomes:

92
 Calculus of Variations (Adapted from Ref. 7)

But, recall that (0) = 0 and (L) = 0, so: 0

So the expression simplifies to:

93
 Calculus of Variations (Adapted from Ref. 7)

Consider the quantity inside the brackets:

We force e to be zero so that:

94
 Calculus of Variations (Adapted from Ref. 7)

Since (x) is arbitrary, this forces us to conclude that

the bracketed quantity:

Note that since e = 0, the ũ and ũ' become u and u'.

This is called the Euler-Lagrange equation and we

will use it to establish the functional for heat
transfer.
95
 Forming the Functional

Recall the governing differential equation for steady

state, one-dimensional conduction:

Area, A, cancels out and we can re-write this as:

where

96
 Forming the Functional

Next, we note that T' replaces u' and which means

that T replaces u in the Euler-Lagrange equation:

So the equation becomes:

97
 Forming the Functional

The heat equation can be rewritten as:

In comparing it with our expression:

We see that:

98
 Forming the Functional

And since:

We can integrate this expression with respect to T'

to get:
0

99
 Forming the Functional

Finally, we arrive at the desired expression:

100
 Finite Element Method Strategy

We will employ the following strategy to solve a

sample problem:

Form the variational statement;

Formulate relations at the element level;
Minimize the integral at the element level;
Assemble equations into a system for solution;
Apply boundary conditions;
Solve for nodal temperatures.

101
 Example 4: One-Dimensional Steady State
Finite Element
Consider, the following configuration:

20 C 5000 W/m 20 C

L = 0.02 m

L = 0.1 m

What is the steady state temperature distribution in

the rod?
102
 Example 4: One-Dimensional Steady State
Finite Element
This is a steady state conduction problem with heat
application/generation and fixed boundary
temperatures;

The governing differential equation is:

0

but for steady state, the right hand side is equal to

zero.
103
 Example 4: One-Dimensional Steady State
Finite Element
First, we form the variational relationship;

From the Euler-Lagrange equation, we have:

And note that:

104
 Example 4: One-Dimensional Steady State
Finite Element
So the Euler-Lagrange equation becomes:

We can rewrite our governing differential equation

in this form:

105
 Example 4: One-Dimensional Steady State
Finite Element
By comparing terms, we see that:

and...

106
 Example 4: One-Dimensional Steady State
Finite Element
We can solve for F by integrating the first
expression:

But, by integrating the second expression (i.e., the

one with the T' term), we get another expression for
F.

107
 Example 4: One-Dimensional Steady State
Finite Element
But both of these expressions must be equal to one
another because they both represent F;

So we conclude that...

And the integral we seek to minimize is:

108
 Example 4: One-Dimensional Steady State
Finite Element (Adapted from Ref. 7)
Now that we have our variational statement, we
need to apply it to our problem;

Most geometries are quite complex and we cannot

apply this expression directly;

So, we discretize the system into finite elements and

apply the expression to each element:

109
 Example 4: One-Dimensional Steady State
Finite Element
Let's consider the following element breakdown;
5000 W/m

20 C (1) (2) (3) (4) (5) 20 C

1 2 3 4 5 6

Note that for finite elements, the node points are at

the boundary of each element;

For finite difference, the nodes were centered in the

element.
110
 Example 4: One-Dimensional Transient Finite
Element (Adapted from Ref. 7)
In this case, we seek to minimize an integral
expression considering the effects of conduction (Ik),
and applied (or internally generated) heating (Iq):

111
 Example 4: One-Dimensional Steady State
Finite Element (Adapted from Ref. 7)
For an individual element, our integral becomes:

where xi and xj represent the bounds of the element,

A(e) is the element cross sectional area and is the
heating per unit length per unit time.

i (e) j
xi xj
x
112
 Example 4: One-Dimensional Steady State
Finite Element (Adapted from Ref. 7)
Let's assume a linear temperature profile across
each element

At the element boundaries, we have:

But we are left with unknown constants, c1 and c2.

113
 Example 4: One-Dimensional Steady State
Finite Element (Adapted from Ref. 7)
But, with two equations and two unknowns, we can
solve for constants c1 and c2.

114
 Example 4: One-Dimensional Steady State
Finite Element (Adapted from Ref. 7)
But, we can express the unknown constants in terms
of known variables from the Ti and Tj expressions to
give us an expression for T(e) in terms of, either,
known or to-be-determined quantities:

115
 Example 4: One-Dimensional Steady State
Finite Element (Adapted from Ref. 7)
We'll also need an expression for the derivative of T
with respect to x;

Differentiating the previous expression yields:

116
 Example 4: One-Dimensional Steady State
Finite Element (Adapted from Ref. 8)
Similarly, we must establish an expression for the
distribution of heating across the element:

As with our temperature calculation, this leads to:

117
 Example 4: One-Dimensional Steady State
Finite Element (Adapted from Ref. 7)
So, starting with the integral expression for an
element, we have:

118
 Example 4: One-Dimensional Steady State
Finite Element (Adapted from Ref. 7)
And substituting our expressions for the element
temperature and gradient, we arrive at:

119
 Example 4: One-Dimensional Steady State
Finite Element (Adapted from Ref. 7)
Before we integrate, we differentiate this previous
expression with respect to Ti :

120
 Example 4: One-Dimensional Steady State
Finite Element (Adapted from Ref. 7)
And we arrive at...

121
 Example 4: One-Dimensional Steady State
Finite Element (Adapted from Ref. 7)
We can now integrate the expression with respect
to x:

122
 Example 4: One-Dimensional Steady State
Finite Element (Adapted from Ref. 7)
This simplifies to:

But, for constant heating across an element:

So the expression becomes:

123
 Example 4: One-Dimensional Steady State
Finite Element (Adapted from Ref. 7)
Remember, the overall integral, I, is a function of all
of the m temperatures in the network:

To minimize the overall integral, we'll need to find

the derivative of each element integral with respect
to every temperature and set the sum equal to zero:

124
 Example 4: One-Dimensional Steady State
Finite Element (Adapted from Ref. 7)
For our linear, one-dimensional elements, though,
I(e) is a function of only two temperatures;

The expression reduces to:

(a) (b)
m-1 m m+1

125
 Example 4: One-Dimensional Steady State
Finite Element (Adapted from Ref. 7)
To simplify matters for our example, let's assume all
elements are the same size, so:

First, for conduction only, the expression reduces to:

(a) (b)
m-1 m m+1

126
 Example 4: One-Dimensional Steady State
Finite Element (Adapted from Ref. 7)
The expression is different when a boundary is
considered;

For our example, both the left and right boundary

temperatures are the same, Tbound :

127
 Example 4: One-Dimensional Steady State
Finite Element (Adapted from Ref. 7)
We can perform a similar operation for heating and
we arrive at:

As a result, half of the heating is applied to one

node and the remaining half to the other.

128
 Example 4: One-Dimensional Steady State
Finite Element (Adapted from Ref. 7)
The system of equations in matrix form is given by:

Not used

Not used

129
 Example 4: One-Dimensional Steady State
Finite Element (Adapted from Ref. 7)
But, before solving, we must apply the boundary
conditions

130
 Example 4: One-Dimensional Steady State
Finite Element (Adapted from Ref. 7)
Moving the terms associated with the boundary
temperatures over to the right hand side yields a
reduced matrix:

131
 Example 4: One-Dimensional Steady State
Finite Element
We see this is of the form:

where...

132
 Example 4: One-Dimensional Steady State
Finite Element
For our sample problem, consider 5 equally sized
elements with 6 equally spaced nodes:

L = 0.1 m
r = 0.01 m

k = 167 W/mC
Tboundary = 20 C
.
qL = 5000 W/m
133
 Example 4: One-Dimensional Steady State
Finite Element
This leads to the following derived quantities:

Dx = 0.02 m
A = 3.1415910-4 m2
.
qLDx = 100 W
kA/Dx = 2.623 W/C
Tb1 = Tb6 = 20 C

134
 Example 4: One-Dimensional Steady State
Finite Element
Filling in numbers:

we get...

135
 Example 4: One-Dimensional Steady State
Finite Element
To solve for temperatures:

136
 Example 4: Finite Element Solution
70

60

50
Temperature (C)

Poor discretization
40 artificially flattens
the temperature
30 profile

20

10

0
0 0.02 0.04 0.06 0.08 0.1

X Location (m)
137
 Example 4: Finite Element Solution Revisited

To fix the problem, adjust discretization*;

5000 W/m

20 C (1) (2) (3) (4) (5) (6) (7) (8) (9) (10) 20 C
1 2 3 4 5 6 7 8 9 10 11

We seek a calculation point somewhere within the

heated region to resolve local temperature
differences;

*Note: Discretization was increased across the entire rod to maintain a

constant element size.
138
 Example 4: Finite Element Solution Revisited
70

60

50
Temperature (C)

Improved
40 discretization allows
temperature
variation to be
30 resolved

20

10

0
0 0.02 0.04 0.06 0.08 0.1
X Location (m) 139
 Example 5: One-Dimensional Steady State
Finite Element
Let's solve the Example 2 geometry using the finite
element method;
100 W
20 C (1) (2) (3) (4) 20 C
1 2 3 4 5

Instead of a distributed heat load (as was the case in

Example 4), all heating is applied at the center (i.e.,
Node 3).

140
 Example 5: One-Dimensional Steady State
Finite Element
The system of equations in matrix form is given by:

Not used

Not used

141
 Example 5: One-Dimensional Steady State
Finite Element
But, before solving, we must apply the boundary
conditions

Not used

Not used

142
 Example 5: One-Dimensional Steady State
Finite Element
Moving the terms associated with the boundary
temperatures over to the right hand side yields the
reduced matrix:

143
 Example 5: One-Dimensional Steady State
Finite Element
Consider 4 equally spaced elements with 5 nodes:

L = 0.1 m
r = 0.01 m

k = 167 W/mC
Tboundary = 20 C
.
Q = 100 W

144
 Example 5: One-Dimensional Steady State
Finite Element
This leads to the following derived quantities:

Dx = 0.025 m
A = 3.1415910-4 m2
.
Q = 100 W
kA/Dx = 2.099 W/C
Tb1 = Tb5 = 20 C

145
 Example 5: One-Dimensional Steady State
Finite Element
Filling in numbers:

we get...

146
 Example 5: One-Dimensional Steady State
Finite Element
To solve for temperatures:

147
 Example 5: Finite Element Solution
70

60

50
Temperature (C)

40

30

20

10

0
0 0.02 0.04 0.06 0.08 0.1
X Location (m) 148
 Example 6: One-Dimensional Transient Finite
Element (Adapted from Ref. 7)
Let's look at the same problem geometry but, this
time, we'll model the transient response.

In this case, we seek to minimize an integral

expression considering the effects of conduction (Ik),
applied heating (Iq) and capacitance (Ic):

149
 Example 6: One-Dimensional Transient Finite
Element (Adapted from Ref. 7)
The individual components are:

where is heating rate and is time.

150
 Example 6: One-Dimensional Transient Finite
Element (Adapted from Ref. 7)
And the overall integral to be minimized is:

We already have expressions for the conduction and

heating so let's focus on the capacitance term.

151
 Example 6: One-Dimensional Transient Finite
Element (Adapted from Ref. 7)
Recall, from Example 4, our expression for
temperature as a function of location in an element:

And because we know the temperatures at the

ends, we solved for constants c1 and c2:

152
 Example 6: One-Dimensional Transient Finite
Element (Adapted from Ref. 7)
So, for an individual element, the integral becomes:

After substituting in expressions for c1 and c2, the

expression inside the integral is a function of, only,
Ti, Tj, xi, xj and x.

153
 Example 6: One-Dimensional Transient Finite
Element (Adapted from Ref. 7)
We differentiate the expression with respect to Ti
and then integrate over x to obtain:

Similarly, we perform the differentiation with

respect to Tj and integrate, we obtain:

154
 Example 6: One-Dimensional Transient Finite
Element (Adapted from Ref. 7)
So, for an element, we see that...

(a) (b)
m-1 m m+1

155
 Example 6: One-Dimensional Transient Finite
Element (Adapted from Ref. 7)
But...

and over a specified time interval, Dq, gives us...

156
 Example 6: One-Dimensional Transient Finite
Element (Adapted from Ref. 7)
The previous expressions form the basis for our
element capacitance matrix:

157
 Example 6: One-Dimensional Transient Finite
Element
But we can assemble the element matrices into a
global capacitance matrix through superposition;

For similarly sized elements, the matrix becomes:

158
 Example 6: One-Dimensional Transient Finite
Element
With the previously derived conduction matrix, we
can now assemble our system of equations:

Becomes...

159
 Example 6: One-Dimensional Transient Finite
Element (Ref. 9)
This equation can be expressed in a more general
form :

so when...
f = 0, we have the Forward difference
f = 1/2, we have Crank-Nicolson
f = 2/3, we have Galerkin
f = 1, we have the Backward difference
160
 Example 6: One-Dimensional Transient Finite
Element (Ref. 9)
Let's use Backward differencing to solve this
problem, i.e., f = 1;

When we substitute into the general equation, we

arrive at:

Our unknowns are the temperatures at time n+1,

shown in blue.

161
 Example 6: One-Dimensional Transient Finite
Element (Ref. 9)
We recognize this equation is of the form:

Inverting [A] and multiplying it with {B} yields:

But we must remember to apply the boundary

conditions as we did before.

162
 Example 6: One-Dimensional Transient Finite
Element
70

60

50
Temperature (C)

40
T1 T2 T3
30 T4 T5

20

10
0 50 100 150 200 250 300
Time (s)
163
 Comparing Finite Difference, Finite Element
with Exact Solution
70
Exact
60
FD
Temperature (C)

50 FE

40

30

20

10

0
0 0.02 0.04 0.06 0.08 0.1
X Location (m) 164
 Conclusion

An overview of numerical methods in heat transfer

has been presented;

The governing differential equation was formulated

from first principles;

The finite difference was developed and

demonstrated through examples;

Finite element was developed and demonstrated

through examples.
165
 Acknowledgements

Dr. Vitali Volovoi of the NESC Statistics Team and Dr.

Jeff Rickman of Lehigh University are acknowledged
for their assistance and guidance with some key
calculations.

Ms. Carol Mosier, Mrs. Melissa Flores, Mr. Richard

Wear, Ms. Ruth Amundsen, Mr. Arturo Avila, Mr.
David Gilmore, Dr. Bruce Drolen, Ms. Robin Beck, Mr.
Hai Nguyen, Mr. John Sharp, Mr. Jack Ercol and Mr.
Hume Peabody are acknowledged for their helpful
review of the draft lesson and for their comments.
166
 References/Credits

1. Van Wylen, G. J. and Sonntag, R. E., Fundamentals of Classical Thermodynamics,

SI Version 2e, New York, John Wiley and Sons, 1978.

2. Holman, J. P., Heat Transfer (Fifth Edition), New York, McGraw Hill Book Company,
1981.

3. Anderson, John D., Jr., Computational Fluid Dynamics, The Basics with
Applications, New York, McGraw-Hill, Inc., 1995.

4. Hornbeck, R., Numerical Methods, Englewood Cliffs, Prentice-Hall, Inc., 1975.

5. McMurchy, R., Thermal Network Modeling Handbook, 14690-H003-R0-00, NASA

Contract 9-10435, January 1972.

167
 References/Credits

6. Jiji, L. M., Heat Conduction, Heidelberg, Springer-Verlag, 2009.

7. Myers, G. E., Analytical Methods in Conduction Heat Transfer, Schenectady,

Genium Publishing Corporation, 1987.

8. Akin, J. E., Finite Element for Undergraduates, Academic Press, 1985.

9. Ramamurty, G., Applied Finite Element Analysis, New Dehli, I.K. International
Publishing Houst Pvt. Ltd., 2009.

Microsoft® Clip Art was used in this presentation.

Wolfram Mathematica® was used for some calculations in this presentation.

168
 For Additional Information

Variational Principles

Liu, G., Exact Variational Principle for 3-D Unsteady Heat Conduction with Second
Sound, Journal of Thermal Science, Vol. 15, No. 4, December 2006.

Vujanovic, B. D. and Jones, S. E., Variational Methods in Nonconservative

Phenomena, Boston, Academic Press, Inc., 1989.

Shih, T. M. (editor), Numerical Properties and Methodologies in Heat Transfer,

Proceedings of the Second National Symposium, Washington, Hemisphere
Publishing Corporation, 1983.

169
 To Contact the Author

Address:
Steven L. Rickman
NASA Engineering and Safety Center (NESC)
NASA - Lyndon B. Johnson Space Center
2101 NASA Parkway
Mail Code: WE
Houston, TX 77058

Phone:
281-483-8867

Email:
steven.l.rickman@nasa.gov

170