UCES
Course: Methods and Analysis in UCES
Applied Area: Fourier heat law, heat diffusion in a thin wire
Model: finite difference model, tridiagonal system
Assessment: thermal properties, radial diffusion, mesh size, finite difference errors
Prerequisites:
Objectives:
General Information: This is the first numerical solution of a boundary value problem. The steady state heat conduction via diffusion and via cooling to the surrounding volume is studied. This type of problem is later generalized to two and three space variables as well as more complicated boundary conditions.
Contact Person: R. E. White, NCSU, white@math.ncsu.edu
Revision Date: 8-3-94
Copyright ©1994, R. E. White. All rights reserved.
In the past sections we have considered differential equations whose solutions were dependent on the time variable. The main physical illustration of this was the heat transfer, but we assumed a uniform temperature distribution in the material. This is rarely the case and one can expect to have a varying temperature with respect to space, and this will result in the diffusion of heat. The math model for this is also a differential equation, but at least one of the independent variables must be space. The simplest boundary value problems have the form
| -Kuxx + Cu = f and boundary conditions u(0), u(1) given. | (1) |
Also, there are other types of boundary conditions which reflect how fast heat passes through the boundary.
The derivation of (1) for steady state one space dimension is based on the empirical law known as
|
| Figure: Heat Diffusion in a Thin Wire |
Consider a thin wire so that there is diffusion in only one direction, x. The f has units of (heat)/(volume time). If there is no time dependence for the temperature, then Fourier's heat law for a thin wire when coupled with Newton's law of cooling to the surrounding region takes the form
If r is the radius of the wire and
r2 ,
f(x) hA 
t
2r c(usur - u) t
t - Kux(x)A t.
K is the constant of proportionality and will depend on the type of material in
the wire. If the temperature is larger on the right side of the cylindrical
cell depicted in the above figure, then ux(x+h) will be positive and
heat will flow into the cell from the right side. A similar argument holds for
the left side except the sign will change. This is an approximation because u
changes continuously. The approximation will get better as h decreases.
Divide both sides of this approximation by hA t and let h go to zero
to get a variation of (1)
| -(Kux )x + 2c/r u = f + 2c/r usur . | (2) |
If K, C and f are constants, then the closed form solution of (1) or (2) is
relatively easy to find. However, if they are more complicated or if we have
diffusion in two and three dimensional space, then closed form solutions are
harder to find. An alternative is the finite difference method which is a way
of converting continuum problems such as (1) or (2) into a finite set of
algebraic equations. It uses numerical derivative approximation for the
second derivative that was discussed at the end of the last section. First,
we break the wire into n equal parts with
The finite difference method or discrete
model approximation to (1) is for
| -K[(ui+1 - ui )/h - (ui - ui-1 )/h]/h + Cui = fi = f(ih). | (3) |
The discrete system (3) may be written in matrix form. For ease of notation
we multiply (3) by h2 , let
Bu1 - Ku2 = h2f1 + Ku0
-Ku1 + Bu2 - Ku3 = h2f2
-Ku2 + Bu3 - Ku4 = h2f3
-Ku3 + Bu4 = h2f4 + Ku5
The matrix form of this is the tridiagonal system
Here with n = 5 we have
The solution can be obtained by either using the tridiagonal (or Thomas)
algorithm, or using a solver that is provided with your computer software.
Let us consider the tridiagonal system Ax = d where A is
.
The following important algorithm will be discussed in detail in the following sections.
(1) = a(1)
(1)= c(1)/(1)
y(1) = d(1)/(1)
for i = 2, n
(i) = a(i)- b(i)*(i-1)
(i) = c(i)/(i)
y(i) = (d(i) - b(i)*y(i-1))/(i)
endloop
x(n) = y(n)
for i = n - 1,1
x(i) = y(i) - (i)*x(i+1)
endloop.
Here we use the Maple procedure linsolve to solve the resulting algebraic
system. The procedure BVP approximates the solution of the thin wire in (2)
where a heat source
Maple Code for Steady State Heat Transfer in a Thin Wire.
Procedure is Defined:
> BVP:=proc(n,tcond,r,c,usur,f)
> with(linalg):
> with(plots):
> C:=(2/r)*c:
> h:= 1/n:
> Rhs:= h*h*(f + C*usur):
> A:=matrix(n-1,n-1,0):
> F:=vector(n-1,0):
> x:=vector(n-1,0):
> x[1]:=h:
> F[1]:= Rhs:
> A[1,1]:=2*tcond + h*h*C:
> A[1,2]:=-tcond:
> for i from 2 to n-2 do
> x[i]:=i*h:
> F[i]:=Rhs:
> A[i,i-1]:= -tcond:
> A[i,i]:= 2*tcond + h*h*C:
> A[i,i+1]:= -tcond:
> od:
> x[n-1]:=(n-1)*h:
> F[n-1]:=Rhs:
> A[n-1,n-2]:= -tcond:
> A[n-1,n-1]:= 2*tcond +h*h*C:
> solution:=linsolve(A,F):
> end:
Input Data:
> n := 10:
> tcond := .001:
> r := .1:
> c := .01:
> usur := 0:
> f := 1:
Call Procedure and Graph Outputs:
The first output uses the above data and the maximum steady state temperature is about 5.
> BVP(10,.001,.1,.01,0,1);
> L:=[[0,0],seq([x[i],solution[i]],i=1..9),[1,0]]:
Next we increase the radius from .1 to .2 and the resulting maximum steady state temperature is about 10.
> BVP(10,.001,.2,.01,0,1);
> LL:=[[0,0],seq([x[i],solution[i]],i=1..9),[1,0]]:
Now we increase the radius to .3 and the maximum steady state temperature is a little above 14.
> BVP(10,.001,.3,.01,0,1);
> LLL:=[[0,0],seq([x[i],solution[i]],i=1..9),[1,0]]:
> plot({L,LL,LLL},x=0..1);
 |
| Figure: Temperature in Thin Wire, r = .1, .2, .3 |
The last calculation keeps the radius at .3, but decreases the c from .01 to .00001 which decreases the amount of heat lost to the surrounding region. This corresponds to the steady state solution that we previously obtained by using the explicit time discretization.
> BVP(10,.001,.3,.00001,0,1);
 |
| Figure: Temperature in a Thin Wire, c = .00001 |
The following Matlab code uses the tridiagonal algorithm to appproximate the
solution of
Matlab Code for a Simple Boundary Value Problem.
function x = bvp(n)
h = 1/n;
for i = 1:n-1
a(i) = 2 + h*h;
b(i) = -1;
c(i) = -1;
d(i) = h*h*f(i*h);
end
trid(n-1,a,b,c,d);
function trid(n,a,b,c,d)
alpha = zeros(n,1);
gamma = zeros(n,1);
y = zeros(n,1);
alpha(1) = a(1);
gamma(1) = c(1)/alpha(1);
y(1) = d(1)/alpha(1);
for i = 2:n
alpha(i) = a(i) - b(i)*gamma(i-1);
gamma(i) = c(i)/alpha(i);
y(i) = (d(i) - b(i)*y(i-1))/alpha(i);
end
x(n) = y(n);
for i = n-1:-1:1
x(i) = y(i) - gamma(i)*x(i+1);
end
In the above model of heat conduction, the thermal conductivity was held constant relative to the space and temperature. If the temperature varies over a large range, the thermal conductivity will show significant changes. Also, the electrical resistivity will vary with temperature, and hence, the heat source , f, will be a function of temperature.
Another important consideration in this model is the possibility of the temperature being a function of the radial space variable, that is, as the radius increases, the temperature is likely to be vary. Hence, heat will also diffuse in this second space direction.
A third consideration is the choice of mesh size, h.
Once the algebraic system has been solved, one wonders how close the numerical
solution of the finite difference method (3), the discrete model, is to the
solution of the differential equation (1), the continuum model. We want to
analyze the discretization error
| K(u((i-1)h) - 2u(ih) + u((i+1)h))/h2 + f(ih) - Cu(ih) = Ku(4)(ci )/12 h2. | (4) |
The finite difference method (3) gives
| -K(ui-1 - 2ui + ui+1 )/h2 + Cui = f(ih). | (5) |
Add equations (4) and (5) and use the definition of Ei to obtain
Or,
Let E = max | Ei | where i ranges from 1 to
Then for all i
2KE/h2 + KM4 /12 h2 .
So, it must be true for the i that gives the maximum, and therefore,
2KE/h2 +
KM4 /12 h2 and
KM4 /12 h2.
We have just proved the next theorem.
Finite Difference Method Error Theorem.
Consider the solution of (1) and the associated finite difference system (3).
If the solution (1) has four continuous derivatives on [0,1], then for
M4/(12C)
h2 .
Example. This example will illustrate the second order convergence of
the finite difference method which was established in the above theorem.
Consider (1) with
| n (h = 1/n) | max. error = E | E/(h2 ) |
| 10 | 17.0267 10-4 | 0.1703 |
| 20 | 4.2140 10-4 | 0.1686 |
| 40 | 1.0865 10-4 | 0.1738 |
| 80 | 0.2598 10-4 | 0.1663 |
| 160 | 0.0585 10-4 | 0.1498 |