AMME2000/BMET2960/BMET9960 - Assignment 2
Due: 11:59 pm Friday 4th June 2021
A team is using drones to capture high resolution images of the terrain after recent volcanic activity
in Iceland. However, the images are distorted due to variations in air pressure and density caused
by the heating from the lava. The team has provided you with some data you can use to obtain a
model of the pressure and density fields in the environment - which will then allow them to apply
corrections to their drone images.
In this assignment you will use the analytical and numerical techniques developed in this course
to obtain a steady-state temperature field in a 2D rectangular slice of the environment using the
Laplace Equation. Using this information, you will then solve the Poisson equation for pressure.
A schematic of the 2D environment is shown in Figure 1.
L = 25m
T(0, z) = T
T(x,0) = T0
T(L, z) = T
T(x,D) = f(x)
Figure 1: Schematic of the 2D domain.
1 Analytic Solution to the Steady-State Temperature Field (20%)
Taking x to be the horizontal position and z to be the depth (positive downwards) for the 2D
domain, the boundary temperature distributions are defined as (Figure 1):
T (0, z) = T (L, z) = T (x, 0) = T0 (1)
T (x,D) = f(x) = T0 + T1 sin(pix/L) + T2 sin(5pix/L) (2)
where L = 25m is the width of the environment in the x-direction, D = 25m is the depth of the
environment in the z-direction, and the temperatures are T0 = 20
◦C, T1 = 380◦C and T2 = 205◦C.
Given that the temperature field is at steady-state we can use the Laplace equation:
= 0 (3)
1. Assuming that the temperature distribution has the form T (x, z) = F (x)G(z), use separation
of variables to show that:
F (x) = A cos(px) +B sin(px) (4)
G(z) = Cepz +De−pz (5)
2. By recognising that T0 appears on all four boundaries, derive the steady-state solution as:
T (x, z) = T0 +
Show your working and justify any assumptions you make. (10%)
3. Explain why Fourier coefficient solutions only exist for n = 1 and n = 5, and hence state the
exact analytic solution for the steady-state temperature field. (5%)
2 Numerical Solution to the Temperature Field (40%)
In this section you will solve the temperature field numerically using two different methods: (i)
Gauss-Seidel and (ii) Successive Over-Relaxation (SOR). We can discretise the spatial domain into
nodes along x and z according to:
xi = (i− 1)∆x where i ∈ [1, nx] (7)
zj = (j − 1)∆z where j ∈ [1, nz] (8)
and denote the temperature at position (xi, zj) by Ti,j . Because we are considering a steady-state
problem, time is not a variable, however we still need to iterate our numerical solution and so will
use the superscript n to index the current solution iteration. In other words, Tni,j represents the
numerical temperature estimate at location (xi, zj) for iteration n.
The Gauss-Seidel (G-S) method computes successive grid point values according to:
Tni+1,j + T
i−1,j + T
i,j+1 + T
Here, provided that you maintain an array Tni,j for the current iteration and a second array T
for the next iteration, you can use this scheme to sequentially compute all the component values
Tn+1i,j . This concept will be illustrated in the lecture.
The Successive Over-Relaxation (SOR) method is a modified form of the G-S method which
computes successive grid point values according to:
Tn+1i,j = T
i,j + 1.85
Tni+1,j + T
i−1,j + T
i,j+1 + T
Using this information, complete the following tasks:
1. Using the information above, explain how you will implement these numerical methods in
MATLAB. In your explanation, you should either include a control flow diagram or steps
with pseudo code. Make sure to outline how you will:
• Implement the boundary conditions
• Initialise the solution domain
• Structure your code to iterate and update the solution
• Check for convergence
2. Solve the steady-state temperature field for the given boundary conditions using both the G-S
and SOR methods. Choose the number of divisions in the domain to be ndiv = 2
n for n = 2 : 9
in both the x and z directions. (Note: this means the number of points is npts = ndiv + 1.)
For each method and each domain discretisation, compute the number of iterations taken to
reach the steady-state solution within a tolerance of = 1× 10−6. Present a figure showing
the number of iterations vs grid size for the two methods. (15%)
3. Using the SOR method, produce three heat-maps showing the solution resolution improve-
ment at 3 different grid sizes tested above (choose a coarse, medium and fine grid to illustrate
the changes). (5%)
Tip: In order to produce a heat-map in MATLAB, use the surf() function with edge
colour set to zero. The syntax for this is: surf(X,Y,T,’EdgeColor’,’none’) where: X is a
matrix of x positions; Y is a matrix of y positions; and T is your 2D temperature solution.
Use the command view(0,-90) to project the surface into the x, y plane, with the colour
gradient representing the temperature in that cell. The command colormap() allows you to
change the colour map, see https://au.mathworks.com/help/matlab/ref/colormap.html
for more options.
4. Compute the L1, L2 and L∞ error norms and the order of the error for each respective norm
for the SOR method (you don’t need to do this for the G-S method). For this error norm
analysis, use the same number of points tested in Part 2.2. The error you compute should be
based on a comparison with the analytic solution obtained in Section 1. Present your results
in a table, as in the lecture slides. State the number of divisions in the domain required to
achieve asymptotic convergence. (10%)
The Lp norm for a two-dimensional solution can be computed using the expression:
|Ti,j − T (xi, zj)|p
Note when computing the norm in MATLAB, we need to supply the norm() function with
a vector input. You can use the following code to compute the p-norms for your solution:
% Create a variable for the difference between solutions:
diff = T numeric(2:end-1,2:end-1) - T analytic(2:end-1,2:end-1);
% "Unwrap" the 2D array to a 1D row vector;
diffVec = reshape(diff',1,numel(diff));
% Compute Lp norms as usual;
L1norm = norm(diffVec,1)*(dx*dz);
L2norm = norm(diffVec,2)*sqrt(dx*dz);
LInfnorm = norm(diffVec,Inf);
3 Poisson Equation for Pressure-Density Coupling (20%)
In this section you will extend the applicability of your SOR solver from Part 2 to the Poisson
Equation. Ultimately, we want to compute the density field, however the density field depends on
both the temperature (which you computed in Part 2) and hydrostatic pressure. The hydrostatic
pressure field can be calculated by solving a Poisson equation:
= φ(x, z) (12)
where ρ is the density and P is the pressure. Given the pressure field from the above equation,
and temperature field from the previous section, the density field can be computed using the Ideal
where Rspecific = 287.058Jkg
−1K−1 is the specific gas constant, and T is the temperature (specified
in Kelvin for this equation).
To solve the Poisson equation, the SOR method may be adapted to include the source term φ(x, z)
Pn+1i,j = P
i,j + 1.85
Pni+1,j + P
i−1,j + P
i,j+1 + P
You can take the boundary conditions of the pressure field to be:
P (x, z = 0) = 101325Pa (15)
∂P (x = 0, z)
∂P (x = L, z)
= 0 (16)
∂P (x, z = D)
= −ρ(x, z = D)g (17)
where gravitational acceleration g = 9.81m/s2. These boundary conditions can be discretised as
Pni,j=1 = 101325 (18)
Pni=1,j = P
Pni=nx,j = P
Pni,j=nz = P
i,j=nz−1 + ρ
Note that the boundary conditions for the density field must be computed using the Ideal Gas
Conceptually, the sequence of steps you will follow is:
i To start your simulation, assume the pressure distribution is P = 101325Pa everywhere.
ii Use the initial pressure field together with the computed temperature field to obtain an initial
guess for the density field.
iii With the initial guess for pressure and density, you can now solve for the new value of pressure
after application of the SOR method:
• Compute the derivatives for the source term.
• Using your modified SOR solver, compute the new pressure field Pn+1i,j .
iv Use the new pressure field to update the density field at the computed points (Equation 13).
v Update the boundary point values of pressure based on the new pressure field.
vi Update the density at the boundary points based on the new pressure at the boundary.
vii Repeat steps (iii) to (vi) until the residuals are below a specified tolerance.
Note that a useful check that your solver is giving sensible results is to consider that for a constant
density field, the total pressure difference from top to bottom of the domain should be equal to
Using this information, complete the following tasks:
1. State the discretisation you have used to calculate the source term φ. (5%)
2. Using the information gained in Section 2, and your own grid convergence study for this
problem, determine an appropriate grid size and residual convergence criteria which gives a
grid converged solution. Present two line plots of density and pressure across the mid-height
of the domain (z=12.5 m) for at least three grid resolutions to support your choice. (10%)
3. Using your solver, compute the density and pressure fields for the domain. Produce heat-
(a) Initial guess for the density field;
(b) Final estimate for the density field (i.e. once pressure field has reached convergence);
(c) Final estimate for the pressure field.
Assignment Submission Information
Your assignment is to be presented as a concise report in PDF format. Note that 10% of the
assignment marks are awarded based on the presentation, clarity and functionality of your MAT-
LAB code. A further 10% of the assignment marks are allocated based on overall structure, clarity
and presentation of the report. An electronic copy of the report should be submitted to
Turnitin by the due date. Late submissions will incur a penalty of 5% per day late.
• Marks will be deducted for handwritten answers and screenshots of equations and/or figures.
See Tutorial 2 for detailed notes on how to insert figures into your assignment.
• The report should not exceed 10 pages; additional pages will not be marked, so aim to be
concise with your findings.
• Any figure/table you include in your report must be numbered and must be referred to and
discussed in the main text of your report.
• You are strongly encouraged to read through the ”Example Assignment” we have made
available on Canvas here. This demonstrates the expectations with respect to assignment
layout, explanations and figures/tables.
• You will submit your MATLAB Code independently from your assignment submission, not as
an appendix. Instructions on how to do this will be included in the Assignment 2 Information
page on Canvas.
• You are also required to submit a script, which should solve Section 2.2 for ndiv = 64 using
the SOR method, and return a heat-map of the temperature solution. This code will be run
by the marker and checked as a part of your MATLAB mark.