CH 08
CH 08
CH 08
1 Introduction
Many important problems are encountered in which heat conduction takes place in a medium
inside which energy is being generated by some process. The process in question may the
passage of electical current, a chemical reaction or a nuclear reaction. Conduction in systems
with internal energy generation can take place under transient conditions, like in many
induction heating applications or under steady state conditions as in nuclear engineering
applications. The basics of heat conduction analysis with internal energy generation are
examined below.
−∇ · q + g(r, t) = 0
For a homogeneous material undergoing no phase changes and introducing the definitional
relationship for the heat flux, the general equation above reduces to
∂T
ρCp = ∇ · k∇T + g(r, t)
∂t
1
ρ is the density, Cp is the specific heat, k is the thermal conductivity, and T is the temper-
ature. While the steady state equation becomes
∇ · k∇T + g(r, t) = 0
2
4 Finite Differences Methods
Finally consider now the problem of steady state conduction with internal heat generation
described by Poisson’s equation
∂2T ∂ 2T g(x, y)
+ = −
∂x2 ∂y 2 k
on (x, y) in the interior R = {(x, y) : a < x < b, c < y < d} of a planar region and subject
to
T (x, y) = TS (x, y)
along the boundary S of the region. As long as g and TS are continuous a unique solution
exists.
As before, a simple approach to the numerical solution of this problem consists in parti-
tioning the intervals [a, b] and [c, d] respectively into N − 1 and M − 1 subintervals with step
sizes ∆x = (b − a)/(N − 1) and ∆y = (d − c)/(M − 1) so that the node or mesh point (xi , yj )
is at xi = a + (i − 1)∆x for i = 1, .., N and yj = c + (j − 1)∆y for j = 1, .., M . Using now
second order accurate centered finite difference approximations for the partial derivatives
the following five point formula for the temperature at node (xi , yj ) is obtained.
1 1 1 1 2 2 g(xi , yj )
− 2
Ti−1,j − 2
Ti+1,j − 2
Ti,j−1 − 2
Ti,j+1 + ( 2 + 2
)Ti,j =
∆x ∆x ∆y ∆y ∆x ∆y k
The boundary conditions imposed along the boundary nodes become
3
constant internal heat generation g is described by the following form of the heat or diffusion
equation
d dT
(k ) + g = 0
dx dx
To implement the finite volume method first subdivide the thickness of the slab into a
collection of N adjoining segments of thickness not necessarily of uniform size (finite volumes)
and place a node inside each volume. Thus, an arbitrary node will be called P , its size is
∆x and the nodes to its left and right, respectively W and E.
Note that two different types of nodes result. While interior nodes are surrounded by
finite volume on both sides, boundary nodes contain material only on one side. Here we
shall concentrate on the derivation of discrete equations for the interior nodes. Those for
the boundary nodes will be discussed later.
The distance between nodes W and P is δxw and that between P and E, δxe . The
locations of the finite volume boundaries corresponding to node P will be denoted by w and
e. Finally, the distance between P and e is called δxe− and that between e and E is δxe+ .
If P is located in the center of the finite volume then δxe− = δxe+ = 12 δxe .
To implement the finite volume method we now integrate the above equation over the
span of the finite volume associated with node P , i.e. from xw to xe
Z e Z e
dT dT dT
d(k ) + gdx = (k )|e − (k )|w + g∆x = 0
w dx w dx dx
Next, approximate the derivatives by piecewise linear profiles to give
TE − TP TP − TW
ke − kw + g∆x = 0
δxe δxw
where the conductivities at the finite volume faces are calculated as the harmonic means of
the values at the neighboring nodes, i.e.
1 − δxe+ /δxe δxe+ /δxe −1
ke = [ + ]
kP kE
and
1 − δxw+ /δxw δxw+ /δxw −1
kw = [ + ]
kW kP
Rearranging one obtains the algebraic equation
aP TP = aE TE + aW TW + b
4
kw
aW =
δxw
aP = aE + aW
and
b = g∆x
One algebraic equation like the above, relating the values of T at three contiguous nodes,
is obtained for each of the N nodes. Adding the this set the discrete equations associated
with the boundary nodes one obtains a consistent set of interlinked simultaneous algebraic
equations which must be solved to give the values of T for all nodal locations..
For the special case of constant thermal properties and finite volumes of uniform size
(δxe = δxw = ∆x, the above is easily rearranged as
TE − 2TP + TW g
+ =0
∆x2 k
which coincides with the FD formula obtained before.
d2 T
− = 20x3
dx2
subject to the boundary conditions T (0) = T (1) = 0.
Now we find an approximate solution using the Galerkin finite element method with three
equal size elements and the following global finite element basis functions:
(
1 − 3x for x ∈ [0, 1/3]
φ1 (x) =
0 elsewhere
3x for x ∈ [0, 1/3]
φ2 (x) = 2 − 3x for x ∈ [1/3, 2/3]
0 elsewhere
5
3x − 1
for x ∈ [1/3, 2/3]
φ3 (x) = 3 − 3x for x ∈ [2/3, 1]
0 elsewhere
(
3x − 2 for x ∈ [2/3, 1]
φ4 (x) =
0 elsewhere
The variational representation of the problem is readily obtained by first multiplying the
given differential equation by a function v and subsequently integrating from x = 0 to x = 1.
After integration by parts, the result is
a(T, v) = (f, v)
where
Z 1 dT dv
a(u, v) = dx
0 dx dx
and
Z 1
(f, v) = 20x3 vdx
0
Now, let us proceed to the implementation of the finite element method. Let the fi-
nite element approximation be uh (x) ≈ T (x). Introduce now for simplicity three elements
e1 , e2 , e3 and four nodes i = 1, 2, 3, 4. The nodal values of u are ui , i = 1, 2, 3, 4 and since u
is specified at x = 0 and x = 1, then u1 = u4 = 0 and the only unknowns are u2 and u3 .
The representation of the trial function u in terms of the nodal values then becomes
4
" #
X u2
uh = ui φi = u2 φ2 + u3 φ3 = [φ2 , φ3 ]
i=1
u3
In the Galerkin method the test function v is given directly in terms of the basis functions.
In the present case, since u1 = u4 = 0, v is simply given by
" #
φ2
vh = φ2 + φ3 = [1, 1]
φ3
and
" #
F2
(f, v) = [1, 1]
F3
6
where
Z 1 Z 2/3
dφ2 2 dφ2 2
k22 = ( ) dx = ( ) dx = 6
0 dx 0 dx
Z 1 Z 2/3
dφ2 dφ3 dφ2 dφ3
k23 = k32 = dx = dx = −3
0 dx dx 1/3 dx dx
Z 1 Z 1
dφ3 2 dφ3 2
k33 = ( ) dx = ( ) dx = 6
0 dx 1/3 dx
Z 1 Z 1/3 Z 2/3 4 20 10
F2 = 20x3 φ2 dx = 20x3 3xdx + 20x3 (2 − 3x)dx = + =
0 0 1/3 81 81 27
and
Z 1 Z 2/3 Z 1 49 131 60
F3 = 20x3 φ3 dx = 20x3 (3x − 1)dx + 20x3 (3 − 3x)dx = + =
0 1/3 2/3 81 81 27
Therefore, the finite element equation is
" #" # " #
10
6 −3 u2 27
Ku = F = = 60
−3 6 u3 27