ELEC0030 Numerical Methods 2019
Coursework Problem 2

Exercise 1 (10%)
Consider a finite difference grid of regular spacing x∆ such that the boundary of the region of
interest lies half way between grid lines. If c is a grid point next to the boundary, (as for example
point 171 in the figure below), and a is one just outside of the boundary, derive an expression for
the value of the function, u(x) at point a to use in the equation for the point c of the grid (for
which a is a neighbour), when the condition that must be satisfied at the boundary is:
( )b
b
u h u u
x ∞
∂
= − −
∂

where h is a constant, b indicates the boundary and u∞ is a given value.
Consider that in this case the function varies monotonically between point c and a (there is
not a maximum or a minimum between them, that is, its derivative doesn't change sign there).
Exercise 2 (75%)
Heat transfer in a solid body can be
represented by the equation:
2u u
t
κ∂ = ∇
∂
(1)
with some appropriate boundary conditions
where u is the temperature and κ is the
thermal diffusivity.
The figure shows the xy-cross section of
a solid material body of a large dimension
along the z axis with a regular finite
differences grid superimposed. The
material parameters and other conditions of
the problem are considered as not varying
in the z-direction so the problem can be
approximated as a 2D problem. On the
bottom surface in the figure, a temperature
T0 is applied a t = 0 and is kept constant at
all times after by a powerful heat source.
The side and top boundaries shown with a
thick black line are insulated and the
insulation can be considered perfect. The
surface portion marked by a dotted line is
open to the air outside where there will be
heat loss by convection1. The distant bulk
average temperature in the air is Ta. Before
the heat source is applied at the bottom,
the temperature of the whole body is uniform and equal to Ta.
Use finite differences and the mesh specified in the figure to find the temperature distribution
as a function of time over the structure.

1 There will also be loss by radiation but we’ll leave that out in this problem.

Finite differences grid
ELEC0030 Numerical Methods 2019
Notes:
The points on the bottom surface have not been numbered since their temperature values are
known.
A surface with a perfect insulation means that no heat flux can cross the surface. Heat flux
density is ˆ uu n
n
∂
−∇ ⋅ = −
∂
, where nˆ is the normal to the surface. Then, perfect insulation means:
. .
0
i b
u
n
∂
=
∂
at all points on the insulated boundaries (i.b.).
Convection is the heat transfer due to flow of a fluid. In this case, the air in contact with the
hot surface. The magnitude of this flow will depend of the temperature difference between the
surface and the bulk of the air, that is:
[( ( , ) ]B
B
uk h u x y u
n ∞
∂
− = −
∂
(2)
where h is the convective heat transfer coefficient, k is the thermal conductivity and u∞ is the
average (bulk) air temperature at some distance from the contact surface.
The differential equation (1) above can be normalised to have a non-dimensional function
u(x,y,t) where the constant κ doesn’t appear explicitly ( 2u u
t
∂
= ∇
∂
). Similarly, consider the
normalised boundary condition (2) where h and k in are also one.
Write a Matlab program that will do all the calculations and draw the plots.
The frontal numbering scheme shown in the figure is convenient for the solution of the matrix
problem. However, in a problem like this, it is probably more convenient to use row and column
numbers (i and j) referring to the grid during the assembly of the matrices.
Time stepping
Although the simplest form to discretise the time derivative would be a forward difference, it
will lead to instabilities. Consider instead using a Crank-Nicolson approach.
From the discretised equation, you will notice that the matrices depend on the grid spacing
and time step. If these are kept constant, the matrices will not change during the time stepping
process and their assembly can be left out of the time stepping loop. However, short steps are
needed at the start. If these are kept constant, it will mean a very large number of iterations and a
long running time. Consider instead using a variable time stepping approach, gradually increasing
the time step through the calculations. For example, the use of a grid spacing d = 0.05 and a
starting time step of dt = 0.002 is suggested, with an increase of 15% per iteration: dt=dt*1.15.
Stopping the iterations
A tolerance of 810− on the absolute error is suggested to stop the iterations, when the error is
defined as the norm of the difference between the solution vector u in two successive time steps:
1n nu u −− .
Constructing the generic equation
Note that the grid has 20 rows and 10 columns of nodes with unknown temperature values.
The nodes are arranged in a regular rectangular pattern so their number (k) can be written in
terms of their row and column numbers (i) and (j). (Warning: Do not confuse these grid row and
column numbers with those of the resultant matrices! – the order the matrices will be 200: 200
rows and 200 columns). Find an expression to convert between the grid row-column notation to
the global numbering easily. Assign the row numbers to the grid of numbered points from the
bottom (row 1: points 1 – 10, and so on and for the columns, from left to right (col. 1 contain the
points 1, 11, 21, …).
ELEC0030 Numerical Methods 2019
Assembly of the matrices:
Write down in you report the general form of the algebraic equation that results from
applying the differential equation (1) to a generic node in the structure, for example the one
corresponding to node 146, that will generate the values for row 146 (of length 200) of the
matrices A and B. Identify all nodes for which this equation applies (internal nodes) by their row
and column numbers in the grid and implement that in your program to fill the matrices.
Consider then all special cases one by one: Nodes next to the Dirichlet boundary (2 – 9), those
next to the side boundaries (11, 21, … 121, 20, 30, 40, …190), those next to the top boundary (192
– 199). For convenience, treat the corner nodes (1, 10, 191, 200 – where two boundary
conditions apply), separately.
List these different cases in your report giving the corresponding modified equation. You will
notice that when deriving the discretised form of the equation for all these cases a matrix
equation can be written of the form:
1n n+ = +Au Bu V
Write a piece of code to fill the corresponding values in each case, separating each segment by
comment lines indicating what is intended.
Notice that the matrices are sparse and have the same sparsity pattern. Sparse storage is then
strongly recommended (top marks will only be available if this route is followed). In this case, you
will have to fill arrays for the row and column indices and for the corresponding matrix values.
(Look for the Matlab documentation for the command “sparse”). The vector V will still be defined
as full and needs to be initialised with: V=zeros(200,1). If you decide to assemble full matrices
instead (but you won’t get full marks), start by pre-allocating space for the matrices A and B as:
A=zeros(200,200) to overwrite later with the values that are not zero.
Once the matrices and the vector are assembled, use the command spy(A) to examine the
sparsity pattern of the matrices. It should be the same for both matrices, verify this and copy this
figure in your report. Does this pattern agree with the expected shape, given the expressions used
to calculate the temperature? Comment on this, referring also to the total number of places in
the matrix and the number of nonzero values.
Run the program with the values:
Time step: ∆t = 0.002 (initial value)
Grid spacing (equal for x and y) d: = 0.05
Initial temperature and bulk temperature of the air: Ta = 10°
Heat source temperature: T0 = 80°.
Output of results and plots
Apart from the sparsity pattern, plot the variation of temperature versus time as the iterations
progress. For this, select a sample of points: For example: 5, 15, 25, … That is, those in column 5,
augmented with the boundary values. That will give a sequence of line plots showing the
temperature profile along that line during the iterations. Label the axes and put a title.
At the end of the iterations, draw a surface plot of the final temperature distribution. For this,
you will have to create a 2D array from the results vector u using the Matlab command “reshape”
and augment that array with extra rows and columns to include the values at the boundaries.
Produce also a contour plot of the final temperature distribution over the cross-section of the
structure.
Exercise 3 (15%)
Consider the same situation described in Exercise 2, but now, we are only interested in a
steady state solution. How would eqn. 1 change in this situation? (Hint: The clue is in the name:
ELEC0030 Numerical Methods 2019
steady state). Write down the discretised equation for a generic node in this case and explain how
all the other equations would change.
Modify a copy of the program written for Exercise 2 to solve this case.
Compare the solution you get now with the plots of the results found at the end of the time
stepping in Exercise 2.


Submit your report through the Moodle page as a pdf file. Write your name in the front page.
The report should be a detailed account of the procedure chosen for these calculations, the
results obtained and the comments or discussion requested. It should include the plots asked and
the programs used for the calculations. You may refer to and cite parts of the code in the text if
that is necessary to explain something but submit the complete program for each exercise as an
appendix to the report.
Do not include Introduction, summary, conclusions or theory sections; just report what you
have done including the equations and necessary derivations and include as appendices the
programs used for the three exercises.
Marks will be allocated to the correctness and completeness of the answers, the presentation
of your report and the quality and clarity of the programs written.