QUESTION 3

Consider a flat hollow cylinder of inner radius π cm and outer radius 2π cm. We will consider a case which is axisymmetric (i.e. does not depend on q) and is shown in Figure 1.

Figure 1: A hollow cylinder (left) modelled as a 2D axisymmetric geometry (right)

The temperature distribution T(r, z,t) across the cylinder may be described using the axisymmetric heat equation given in Eqn (3). The cylinder has the following boundary condition:

Q(r, z,t) = sin(2r)   Eqn (4)

and initial condition:

Q(r, z, 0) = sin(2r)   Eqn (5)

Q3a In your PDF, outline the algorithm needed to solve Eqn 3 using the BTCS scheme from time = 0 to time = tn, with a constant timestep, Dt (assume that tn = NDt). Include the major steps required to set the problem up (don’t include minor steps like setting grid spacing, etc.) This algorithm must be in the form of pseudo code or a flow chart. An example of pseudo code for solving an ODE using the Euler method is shown in the snippet below.

Set up parameters of the problem

Determine stepsize and number of steps

for each time step

 Solve the Euler method

 Store calculated values in vector

end

NOTE: All BCs in this problem are Dirichlet conditions (i.e. QBC is known in advance). THUS, it is possible to remove the Boundary points from the set of unknowns, which is what is recommended below. If you want to keep them as unknowns, then the size of the grid will be Nr*Nz instead of (Nr-2)*(Nz-2) as specified in Q3b.

Q3b Write a Matlab function that sets up the Laplacian matrix L and the boundary condition vector BC. Your function header must be written as

function [L,BC] = MatSet(Nr,Nz,rg,zg,alpha,BCf)

 where

• Nr and Nz are the number of points in the r- and z-directions, respectively

• rg and zg are column vectors containing the r- and z-values of the grid points respectively.

• a is the thermal diffusion coefficient

• BCf is a function describing the boundary condition as a function of r,z and t (see Eqn 4)

• L is the ((Nr-2)*(Nz-2) x (Nr-2)*(Nz-2)) Laplacian matrix

• BC is a ((Nr-2)*(Nz-2) x 1) column vector on the RHS that accounts for the boundary conditions

NOTE: You can define an anonymous function that converts a 2D node number (j,k) = (jdr, kdz) to a global node number for the unknowns. Your function will look like

nn=@(j,k,Nr) {some function of j, k, Nr};

(See Workshops 25-27 for examples of what this looks like). When setting up the matrix problem the function nn (or at least a similar process) can be used to reliably determine the row and column number of each entry in the Laplacian matrix and the row number in the BC matrix.

Q3c Modify Lab06_Q3c.m to solve the 2D axisymmetric heat equation using the BTCS scheme (i.e. Implement your algorithm from Q3a) together with the function you have written in Q3b. Use the following values of the parameters

 • a = 0.05,

• Nr = Nz =N (use 4 values 11, 21, 31 and 41),

• time step Dt = 1 s

• solve from t = 0 to t = 50 s for each value of N

Use the commands‘tic’ and ‘toc’ to determine how much time is spent setting the matrix problem up, (i.e. the time spent assembling the Laplacian matrix and BC vector) and separately, how much time is spent undertaking the time-stepping solution. (NOTE: The time taken for plotting should not be included as part of the time taken for time-stepping).

Print out the time required for setting up each grid and solving each solution.

Write a few sentences comparing the times and explain why you think you get the results you get.

Using the subplot command, plot contours of your solution at time=50 including the boundary values for each (Nr, Nz) in a 2x2 grid. You should have one figure window displaying the contours as:

Note that you will have to reshape your 1D vector of unknowns into a 2D matrix and embed that into a slightly bigger 2D matrix so that you can add boundary conditions.

Repeat this question for a time step of time step Δt = 50 s.

Write a few sentences to the command window that compares the results for the 2 different timesteps. Are you surprised at your results for a single timestep of 50s compared to those for multiple steps with a 1 s timestep?

data.all.m